1 /*$Id: aij.c,v 1.372 2001/05/29 21:02:54 bsmith Exp balay $*/ 2 /* 3 Defines the basic matrix operations for the AIJ (compressed row) 4 matrix storage format. 5 */ 6 7 #include "petscsys.h" /*I "petscmat.h" I*/ 8 #include "src/mat/impls/aij/seq/aij.h" 9 #include "src/vec/vecimpl.h" 10 #include "src/inline/spops.h" 11 #include "src/inline/dot.h" 12 #include "petscbt.h" 13 14 15 EXTERN int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 16 17 #undef __FUNCT__ 18 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 19 int MatGetRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *m,int **ia,int **ja,PetscTruth *done) 20 { 21 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 22 int ierr,i,ishift; 23 24 PetscFunctionBegin; 25 *m = A->m; 26 if (!ia) PetscFunctionReturn(0); 27 ishift = a->indexshift; 28 if (symmetric && !A->structurally_symmetric) { 29 ierr = MatToSymmetricIJ_SeqAIJ(A->m,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 30 } else if (oshift == 0 && ishift == -1) { 31 int nz = a->i[A->m] - 1; 32 /* malloc space and subtract 1 from i and j indices */ 33 ierr = PetscMalloc((A->m+1)*sizeof(int),ia);CHKERRQ(ierr); 34 ierr = PetscMalloc((nz+1)*sizeof(int),ja);CHKERRQ(ierr); 35 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] - 1; 36 for (i=0; i<A->m+1; i++) (*ia)[i] = a->i[i] - 1; 37 } else if (oshift == 1 && ishift == 0) { 38 int nz = a->i[A->m]; 39 /* malloc space and add 1 to i and j indices */ 40 ierr = PetscMalloc((A->m+1)*sizeof(int),ia);CHKERRQ(ierr); 41 ierr = PetscMalloc((nz+1)*sizeof(int),ja);CHKERRQ(ierr); 42 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 43 for (i=0; i<A->m+1; i++) (*ia)[i] = a->i[i] + 1; 44 } else { 45 *ia = a->i; *ja = a->j; 46 } 47 PetscFunctionReturn(0); 48 } 49 50 #undef __FUNCT__ 51 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 52 int MatRestoreRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia,int **ja,PetscTruth *done) 53 { 54 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 55 int ishift = a->indexshift,ierr; 56 57 PetscFunctionBegin; 58 if (!ia) PetscFunctionReturn(0); 59 if ((symmetric && !A->structurally_symmetric) || (oshift == 0 && ishift == -1) || (oshift == 1 && ishift == 0)) { 60 ierr = PetscFree(*ia);CHKERRQ(ierr); 61 ierr = PetscFree(*ja);CHKERRQ(ierr); 62 } 63 PetscFunctionReturn(0); 64 } 65 66 #undef __FUNCT__ 67 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 68 int MatGetColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja,PetscTruth *done) 69 { 70 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 71 int ierr,i,ishift = a->indexshift,*collengths,*cia,*cja,n = A->n,m = A->m; 72 int nz = a->i[m]+ishift,row,*jj,mr,col; 73 74 PetscFunctionBegin; 75 *nn = A->n; 76 if (!ia) PetscFunctionReturn(0); 77 if (symmetric) { 78 ierr = MatToSymmetricIJ_SeqAIJ(A->m,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 79 } else { 80 ierr = PetscMalloc((n+1)*sizeof(int),&collengths);CHKERRQ(ierr); 81 ierr = PetscMemzero(collengths,n*sizeof(int));CHKERRQ(ierr); 82 ierr = PetscMalloc((n+1)*sizeof(int),&cia);CHKERRQ(ierr); 83 ierr = PetscMalloc((nz+1)*sizeof(int),&cja);CHKERRQ(ierr); 84 jj = a->j; 85 for (i=0; i<nz; i++) { 86 collengths[jj[i] + ishift]++; 87 } 88 cia[0] = oshift; 89 for (i=0; i<n; i++) { 90 cia[i+1] = cia[i] + collengths[i]; 91 } 92 ierr = PetscMemzero(collengths,n*sizeof(int));CHKERRQ(ierr); 93 jj = a->j; 94 for (row=0; row<m; row++) { 95 mr = a->i[row+1] - a->i[row]; 96 for (i=0; i<mr; i++) { 97 col = *jj++ + ishift; 98 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 99 } 100 } 101 ierr = PetscFree(collengths);CHKERRQ(ierr); 102 *ia = cia; *ja = cja; 103 } 104 PetscFunctionReturn(0); 105 } 106 107 #undef __FUNCT__ 108 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 109 int MatRestoreColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia,int **ja,PetscTruth *done) 110 { 111 int ierr; 112 113 PetscFunctionBegin; 114 if (!ia) PetscFunctionReturn(0); 115 116 ierr = PetscFree(*ia);CHKERRQ(ierr); 117 ierr = PetscFree(*ja);CHKERRQ(ierr); 118 119 PetscFunctionReturn(0); 120 } 121 122 #define CHUNKSIZE 15 123 124 #undef __FUNCT__ 125 #define __FUNCT__ "MatSetValues_SeqAIJ" 126 int MatSetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 127 { 128 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 129 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted = a->sorted; 130 int *imax = a->imax,*ai = a->i,*ailen = a->ilen; 131 int *aj = a->j,nonew = a->nonew,shift = a->indexshift,ierr; 132 Scalar *ap,value,*aa = a->a; 133 PetscTruth ignorezeroentries = ((a->ignorezeroentries && is == ADD_VALUES) ? PETSC_TRUE:PETSC_FALSE); 134 PetscTruth roworiented = a->roworiented; 135 136 PetscFunctionBegin; 137 for (k=0; k<m; k++) { /* loop over added rows */ 138 row = im[k]; 139 if (row < 0) continue; 140 #if defined(PETSC_USE_BOPT_g) 141 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %d max %d",row,A->m); 142 #endif 143 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 144 rmax = imax[row]; nrow = ailen[row]; 145 low = 0; 146 for (l=0; l<n; l++) { /* loop over added columns */ 147 if (in[l] < 0) continue; 148 #if defined(PETSC_USE_BOPT_g) 149 if (in[l] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %d max %d",in[l],A->n); 150 #endif 151 col = in[l] - shift; 152 if (roworiented) { 153 value = v[l + k*n]; 154 } else { 155 value = v[k + l*m]; 156 } 157 if (value == 0.0 && ignorezeroentries) continue; 158 159 if (!sorted) low = 0; high = nrow; 160 while (high-low > 5) { 161 t = (low+high)/2; 162 if (rp[t] > col) high = t; 163 else low = t; 164 } 165 for (i=low; i<high; i++) { 166 if (rp[i] > col) break; 167 if (rp[i] == col) { 168 if (is == ADD_VALUES) ap[i] += value; 169 else ap[i] = value; 170 goto noinsert; 171 } 172 } 173 if (nonew == 1) goto noinsert; 174 else if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%d,%d) in the matrix",row,col); 175 if (nrow >= rmax) { 176 /* there is no extra room in row, therefore enlarge */ 177 int new_nz = ai[A->m] + CHUNKSIZE,len,*new_i,*new_j; 178 Scalar *new_a; 179 180 if (nonew == -2) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%d,%d) in the matrix requiring new malloc()",row,col); 181 182 /* malloc new storage space */ 183 len = new_nz*(sizeof(int)+sizeof(Scalar))+(A->m+1)*sizeof(int); 184 ierr = PetscMalloc(len,&new_a);CHKERRQ(ierr); 185 new_j = (int*)(new_a + new_nz); 186 new_i = new_j + new_nz; 187 188 /* copy over old data into new slots */ 189 for (ii=0; ii<row+1; ii++) {new_i[ii] = ai[ii];} 190 for (ii=row+1; ii<A->m+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 191 ierr = PetscMemcpy(new_j,aj,(ai[row]+nrow+shift)*sizeof(int));CHKERRQ(ierr); 192 len = (new_nz - CHUNKSIZE - ai[row] - nrow - shift); 193 ierr = PetscMemcpy(new_j+ai[row]+shift+nrow+CHUNKSIZE,aj+ai[row]+shift+nrow,len*sizeof(int));CHKERRQ(ierr); 194 ierr = PetscMemcpy(new_a,aa,(ai[row]+nrow+shift)*sizeof(Scalar));CHKERRQ(ierr); 195 ierr = PetscMemcpy(new_a+ai[row]+shift+nrow+CHUNKSIZE,aa+ai[row]+shift+nrow,len*sizeof(Scalar));CHKERRQ(ierr); 196 /* free up old matrix storage */ 197 ierr = PetscFree(a->a);CHKERRQ(ierr); 198 if (!a->singlemalloc) { 199 ierr = PetscFree(a->i);CHKERRQ(ierr); 200 ierr = PetscFree(a->j);CHKERRQ(ierr); 201 } 202 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 203 a->singlemalloc = PETSC_TRUE; 204 205 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 206 rmax = imax[row] = imax[row] + CHUNKSIZE; 207 PetscLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + sizeof(Scalar))); 208 a->maxnz += CHUNKSIZE; 209 a->reallocs++; 210 } 211 N = nrow++ - 1; a->nz++; 212 /* shift up all the later entries in this row */ 213 for (ii=N; ii>=i; ii--) { 214 rp[ii+1] = rp[ii]; 215 ap[ii+1] = ap[ii]; 216 } 217 rp[i] = col; 218 ap[i] = value; 219 noinsert:; 220 low = i + 1; 221 } 222 ailen[row] = nrow; 223 } 224 PetscFunctionReturn(0); 225 } 226 227 #undef __FUNCT__ 228 #define __FUNCT__ "MatGetValues_SeqAIJ" 229 int MatGetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 230 { 231 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 232 int *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 233 int *ai = a->i,*ailen = a->ilen,shift = a->indexshift; 234 Scalar *ap,*aa = a->a,zero = 0.0; 235 236 PetscFunctionBegin; 237 for (k=0; k<m; k++) { /* loop over rows */ 238 row = im[k]; 239 if (row < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %d",row); 240 if (row >= A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: %d",row); 241 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 242 nrow = ailen[row]; 243 for (l=0; l<n; l++) { /* loop over columns */ 244 if (in[l] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %d",in[l]); 245 if (in[l] >= A->n) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: %d",in[l]); 246 col = in[l] - shift; 247 high = nrow; low = 0; /* assume unsorted */ 248 while (high-low > 5) { 249 t = (low+high)/2; 250 if (rp[t] > col) high = t; 251 else low = t; 252 } 253 for (i=low; i<high; i++) { 254 if (rp[i] > col) break; 255 if (rp[i] == col) { 256 *v++ = ap[i]; 257 goto finished; 258 } 259 } 260 *v++ = zero; 261 finished:; 262 } 263 } 264 PetscFunctionReturn(0); 265 } 266 267 268 #undef __FUNCT__ 269 #define __FUNCT__ "MatView_SeqAIJ_Binary" 270 int MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 271 { 272 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 273 int i,fd,*col_lens,ierr; 274 275 PetscFunctionBegin; 276 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 277 ierr = PetscMalloc((4+A->m)*sizeof(int),&col_lens);CHKERRQ(ierr); 278 col_lens[0] = MAT_COOKIE; 279 col_lens[1] = A->m; 280 col_lens[2] = A->n; 281 col_lens[3] = a->nz; 282 283 /* store lengths of each row and write (including header) to file */ 284 for (i=0; i<A->m; i++) { 285 col_lens[4+i] = a->i[i+1] - a->i[i]; 286 } 287 ierr = PetscBinaryWrite(fd,col_lens,4+A->m,PETSC_INT,1);CHKERRQ(ierr); 288 ierr = PetscFree(col_lens);CHKERRQ(ierr); 289 290 /* store column indices (zero start index) */ 291 if (a->indexshift) { 292 for (i=0; i<a->nz; i++) a->j[i]--; 293 } 294 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,0);CHKERRQ(ierr); 295 if (a->indexshift) { 296 for (i=0; i<a->nz; i++) a->j[i]++; 297 } 298 299 /* store nonzero values */ 300 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,0);CHKERRQ(ierr); 301 PetscFunctionReturn(0); 302 } 303 304 #undef __FUNCT__ 305 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 306 int MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 307 { 308 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 309 int ierr,i,j,m = A->m,shift = a->indexshift; 310 char *name; 311 PetscViewerFormat format; 312 313 PetscFunctionBegin; 314 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 315 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 316 if (format == PETSC_VIEWER_ASCII_INFO_LONG || format == PETSC_VIEWER_ASCII_INFO) { 317 if (a->inode.size) { 318 ierr = PetscViewerASCIIPrintf(viewer,"using I-node routines: found %d nodes, limit used is %d\n",a->inode.node_count,a->inode.limit);CHKERRQ(ierr); 319 } else { 320 ierr = PetscViewerASCIIPrintf(viewer,"not using I-node routines\n");CHKERRQ(ierr); 321 } 322 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 323 int nofinalvalue = 0; 324 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->n-!shift)) { 325 nofinalvalue = 1; 326 } 327 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 328 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %d %d \n",m,A->n);CHKERRQ(ierr); 329 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %d \n",a->nz);CHKERRQ(ierr); 330 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%d,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 331 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 332 333 for (i=0; i<m; i++) { 334 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 335 #if defined(PETSC_USE_COMPLEX) 336 ierr = PetscViewerASCIIPrintf(viewer,"%d %d %18.16e + %18.16ei \n",i+1,a->j[j]+!shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 337 #else 338 ierr = PetscViewerASCIIPrintf(viewer,"%d %d %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 339 #endif 340 } 341 } 342 if (nofinalvalue) { 343 ierr = PetscViewerASCIIPrintf(viewer,"%d %d %18.16e\n",m,A->n,0.0);CHKERRQ(ierr); 344 } 345 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 346 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 347 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 348 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 349 for (i=0; i<m; i++) { 350 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i);CHKERRQ(ierr); 351 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 352 #if defined(PETSC_USE_COMPLEX) 353 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 354 ierr = PetscViewerASCIIPrintf(viewer," %d %g + %g i",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 355 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 356 ierr = PetscViewerASCIIPrintf(viewer," %d %g - %g i",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 357 } else if (PetscRealPart(a->a[j]) != 0.0) { 358 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 359 } 360 #else 361 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %d %g ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 362 #endif 363 } 364 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 365 } 366 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 367 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 368 int nzd=0,fshift=1,*sptr; 369 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 370 ierr = PetscMalloc((m+1)*sizeof(int),&sptr);CHKERRQ(ierr); 371 for (i=0; i<m; i++) { 372 sptr[i] = nzd+1; 373 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 374 if (a->j[j] >= i) { 375 #if defined(PETSC_USE_COMPLEX) 376 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 377 #else 378 if (a->a[j] != 0.0) nzd++; 379 #endif 380 } 381 } 382 } 383 sptr[m] = nzd+1; 384 ierr = PetscViewerASCIIPrintf(viewer," %d %d\n\n",m,nzd);CHKERRQ(ierr); 385 for (i=0; i<m+1; i+=6) { 386 if (i+4<m) {ierr = PetscViewerASCIIPrintf(viewer," %d %d %d %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5]);CHKERRQ(ierr);} 387 else if (i+3<m) {ierr = PetscViewerASCIIPrintf(viewer," %d %d %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4]);CHKERRQ(ierr);} 388 else if (i+2<m) {ierr = PetscViewerASCIIPrintf(viewer," %d %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]);CHKERRQ(ierr);} 389 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 390 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %d %d\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 391 else {ierr = PetscViewerASCIIPrintf(viewer," %d\n",sptr[i]);CHKERRQ(ierr);} 392 } 393 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 394 ierr = PetscFree(sptr);CHKERRQ(ierr); 395 for (i=0; i<m; i++) { 396 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 397 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %d ",a->j[j]+fshift);CHKERRQ(ierr);} 398 } 399 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 400 } 401 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 402 for (i=0; i<m; i++) { 403 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 404 if (a->j[j] >= i) { 405 #if defined(PETSC_USE_COMPLEX) 406 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 407 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 408 } 409 #else 410 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 411 #endif 412 } 413 } 414 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 415 } 416 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 417 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 418 int cnt = 0,jcnt; 419 Scalar value; 420 421 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 422 for (i=0; i<m; i++) { 423 jcnt = 0; 424 for (j=0; j<A->n; j++) { 425 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 426 value = a->a[cnt++]; 427 jcnt++; 428 } else { 429 value = 0.0; 430 } 431 #if defined(PETSC_USE_COMPLEX) 432 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 433 #else 434 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 435 #endif 436 } 437 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 438 } 439 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 440 } else { 441 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 442 for (i=0; i<m; i++) { 443 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i);CHKERRQ(ierr); 444 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 445 #if defined(PETSC_USE_COMPLEX) 446 if (PetscImaginaryPart(a->a[j]) > 0.0) { 447 ierr = PetscViewerASCIIPrintf(viewer," %d %g + %g i",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 448 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 449 ierr = PetscViewerASCIIPrintf(viewer," %d %g - %g i",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 450 } else { 451 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 452 } 453 #else 454 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 455 #endif 456 } 457 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 458 } 459 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 460 } 461 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 462 PetscFunctionReturn(0); 463 } 464 465 #undef __FUNCT__ 466 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 467 int MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 468 { 469 Mat A = (Mat) Aa; 470 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 471 int ierr,i,j,m = A->m,shift = a->indexshift,color; 472 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 473 PetscViewer viewer; 474 PetscViewerFormat format; 475 476 PetscFunctionBegin; 477 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 478 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 479 480 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 481 /* loop over matrix elements drawing boxes */ 482 483 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 484 /* Blue for negative, Cyan for zero and Red for positive */ 485 color = PETSC_DRAW_BLUE; 486 for (i=0; i<m; i++) { 487 y_l = m - i - 1.0; y_r = y_l + 1.0; 488 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 489 x_l = a->j[j] + shift; x_r = x_l + 1.0; 490 #if defined(PETSC_USE_COMPLEX) 491 if (PetscRealPart(a->a[j]) >= 0.) continue; 492 #else 493 if (a->a[j] >= 0.) continue; 494 #endif 495 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 496 } 497 } 498 color = PETSC_DRAW_CYAN; 499 for (i=0; i<m; i++) { 500 y_l = m - i - 1.0; y_r = y_l + 1.0; 501 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 502 x_l = a->j[j] + shift; x_r = x_l + 1.0; 503 if (a->a[j] != 0.) continue; 504 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 505 } 506 } 507 color = PETSC_DRAW_RED; 508 for (i=0; i<m; i++) { 509 y_l = m - i - 1.0; y_r = y_l + 1.0; 510 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 511 x_l = a->j[j] + shift; x_r = x_l + 1.0; 512 #if defined(PETSC_USE_COMPLEX) 513 if (PetscRealPart(a->a[j]) <= 0.) continue; 514 #else 515 if (a->a[j] <= 0.) continue; 516 #endif 517 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 518 } 519 } 520 } else { 521 /* use contour shading to indicate magnitude of values */ 522 /* first determine max of all nonzero values */ 523 int nz = a->nz,count; 524 PetscDraw popup; 525 PetscReal scale; 526 527 for (i=0; i<nz; i++) { 528 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 529 } 530 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 531 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 532 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 533 count = 0; 534 for (i=0; i<m; i++) { 535 y_l = m - i - 1.0; y_r = y_l + 1.0; 536 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 537 x_l = a->j[j] + shift; x_r = x_l + 1.0; 538 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(a->a[count])); 539 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 540 count++; 541 } 542 } 543 } 544 PetscFunctionReturn(0); 545 } 546 547 #undef __FUNCT__ 548 #define __FUNCT__ "MatView_SeqAIJ_Draw" 549 int MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 550 { 551 int ierr; 552 PetscDraw draw; 553 PetscReal xr,yr,xl,yl,h,w; 554 PetscTruth isnull; 555 556 PetscFunctionBegin; 557 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 558 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 559 if (isnull) PetscFunctionReturn(0); 560 561 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 562 xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0; 563 xr += w; yr += h; xl = -w; yl = -h; 564 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 565 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 566 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 567 PetscFunctionReturn(0); 568 } 569 570 #undef __FUNCT__ 571 #define __FUNCT__ "MatView_SeqAIJ" 572 int MatView_SeqAIJ(Mat A,PetscViewer viewer) 573 { 574 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 575 int ierr; 576 PetscTruth issocket,isascii,isbinary,isdraw; 577 578 PetscFunctionBegin; 579 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 580 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 581 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 582 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 583 if (issocket) { 584 if (a->indexshift) { 585 SETERRQ(1,"Can only socket send sparse matrix with 0 based indexing"); 586 } 587 ierr = PetscViewerSocketPutSparse_Private(viewer,A->m,A->n,a->nz,a->a,a->i,a->j);CHKERRQ(ierr); 588 } else if (isascii) { 589 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 590 } else if (isbinary) { 591 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 592 } else if (isdraw) { 593 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 594 } else { 595 SETERRQ1(1,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 596 } 597 PetscFunctionReturn(0); 598 } 599 600 EXTERN int Mat_AIJ_CheckInode(Mat); 601 #undef __FUNCT__ 602 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 603 int MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 604 { 605 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 606 int fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax,ierr; 607 int m = A->m,*ip,N,*ailen = a->ilen,shift = a->indexshift,rmax = 0; 608 Scalar *aa = a->a,*ap; 609 610 PetscFunctionBegin; 611 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 612 613 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 614 for (i=1; i<m; i++) { 615 /* move each row back by the amount of empty slots (fshift) before it*/ 616 fshift += imax[i-1] - ailen[i-1]; 617 rmax = PetscMax(rmax,ailen[i]); 618 if (fshift) { 619 ip = aj + ai[i] + shift; 620 ap = aa + ai[i] + shift; 621 N = ailen[i]; 622 for (j=0; j<N; j++) { 623 ip[j-fshift] = ip[j]; 624 ap[j-fshift] = ap[j]; 625 } 626 } 627 ai[i] = ai[i-1] + ailen[i-1]; 628 } 629 if (m) { 630 fshift += imax[m-1] - ailen[m-1]; 631 ai[m] = ai[m-1] + ailen[m-1]; 632 } 633 /* reset ilen and imax for each row */ 634 for (i=0; i<m; i++) { 635 ailen[i] = imax[i] = ai[i+1] - ai[i]; 636 } 637 a->nz = ai[m] + shift; 638 639 /* diagonals may have moved, so kill the diagonal pointers */ 640 if (fshift && a->diag) { 641 ierr = PetscFree(a->diag);CHKERRQ(ierr); 642 PetscLogObjectMemory(A,-(m+1)*sizeof(int)); 643 a->diag = 0; 644 } 645 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %d X %d; storage space: %d unneeded,%d used\n",m,A->n,fshift,a->nz); 646 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues() is %d\n",a->reallocs); 647 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Most nonzeros in any row is %d\n",rmax); 648 a->reallocs = 0; 649 A->info.nz_unneeded = (double)fshift; 650 a->rmax = rmax; 651 652 /* check out for identical nodes. If found, use inode functions */ 653 ierr = Mat_AIJ_CheckInode(A);CHKERRQ(ierr); 654 PetscFunctionReturn(0); 655 } 656 657 #undef __FUNCT__ 658 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 659 int MatZeroEntries_SeqAIJ(Mat A) 660 { 661 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 662 int ierr; 663 664 PetscFunctionBegin; 665 ierr = PetscMemzero(a->a,(a->i[A->m]+a->indexshift)*sizeof(Scalar));CHKERRQ(ierr); 666 PetscFunctionReturn(0); 667 } 668 669 #undef __FUNCT__ 670 #define __FUNCT__ "MatDestroy_SeqAIJ" 671 int MatDestroy_SeqAIJ(Mat A) 672 { 673 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 674 int ierr; 675 676 PetscFunctionBegin; 677 #if defined(PETSC_USE_LOG) 678 PetscLogObjectState((PetscObject)A,"Rows=%d, Cols=%d, NZ=%d",A->m,A->n,a->nz); 679 #endif 680 if (a->freedata) { 681 ierr = PetscFree(a->a);CHKERRQ(ierr); 682 if (!a->singlemalloc) { 683 ierr = PetscFree(a->i);CHKERRQ(ierr); 684 ierr = PetscFree(a->j);CHKERRQ(ierr); 685 } 686 } 687 if (a->row) { 688 ierr = ISDestroy(a->row);CHKERRQ(ierr); 689 } 690 if (a->col) { 691 ierr = ISDestroy(a->col);CHKERRQ(ierr); 692 } 693 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 694 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 695 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 696 if (a->idiag) {ierr = PetscFree(a->idiag);CHKERRQ(ierr);} 697 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 698 if (a->inode.size) {ierr = PetscFree(a->inode.size);CHKERRQ(ierr);} 699 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 700 if (a->saved_values) {ierr = PetscFree(a->saved_values);CHKERRQ(ierr);} 701 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 702 ierr = PetscFree(a);CHKERRQ(ierr); 703 PetscFunctionReturn(0); 704 } 705 706 #undef __FUNCT__ 707 #define __FUNCT__ "MatCompress_SeqAIJ" 708 int MatCompress_SeqAIJ(Mat A) 709 { 710 PetscFunctionBegin; 711 PetscFunctionReturn(0); 712 } 713 714 #undef __FUNCT__ 715 #define __FUNCT__ "MatSetOption_SeqAIJ" 716 int MatSetOption_SeqAIJ(Mat A,MatOption op) 717 { 718 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 719 720 PetscFunctionBegin; 721 if (op == MAT_ROW_ORIENTED) a->roworiented = PETSC_TRUE; 722 else if (op == MAT_KEEP_ZEROED_ROWS) a->keepzeroedrows = PETSC_TRUE; 723 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = PETSC_FALSE; 724 else if (op == MAT_COLUMNS_SORTED) a->sorted = PETSC_TRUE; 725 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = PETSC_FALSE; 726 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 727 else if (op == MAT_NEW_NONZERO_LOCATION_ERR) a->nonew = -1; 728 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERR) a->nonew = -2; 729 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 730 else if (op == MAT_IGNORE_ZERO_ENTRIES) a->ignorezeroentries = PETSC_TRUE; 731 else if (op == MAT_USE_INODES) a->inode.use = PETSC_TRUE; 732 else if (op == MAT_DO_NOT_USE_INODES) a->inode.use = PETSC_FALSE; 733 else if (op == MAT_ROWS_SORTED || 734 op == MAT_ROWS_UNSORTED || 735 op == MAT_YES_NEW_DIAGONALS || 736 op == MAT_IGNORE_OFF_PROC_ENTRIES|| 737 op == MAT_USE_HASH_TABLE) 738 PetscLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 739 else if (op == MAT_NO_NEW_DIAGONALS) { 740 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 741 } else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 742 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 743 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 744 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 745 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 746 else SETERRQ(PETSC_ERR_SUP,"unknown option"); 747 PetscFunctionReturn(0); 748 } 749 750 #undef __FUNCT__ 751 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 752 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 753 { 754 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 755 int i,j,n,shift = a->indexshift,ierr; 756 Scalar *x,zero = 0.0; 757 758 PetscFunctionBegin; 759 ierr = VecSet(&zero,v);CHKERRQ(ierr); 760 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 761 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 762 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 763 for (i=0; i<A->m; i++) { 764 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 765 if (a->j[j]+shift == i) { 766 x[i] = a->a[j]; 767 break; 768 } 769 } 770 } 771 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 772 PetscFunctionReturn(0); 773 } 774 775 #undef __FUNCT__ 776 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 777 int MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 778 { 779 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 780 Scalar *x,*y,*v,alpha,zero = 0.0; 781 int ierr,m = A->m,n,i,*idx,shift = a->indexshift; 782 783 PetscFunctionBegin; 784 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 785 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 786 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 787 y = y + shift; /* shift for Fortran start by 1 indexing */ 788 for (i=0; i<m; i++) { 789 idx = a->j + a->i[i] + shift; 790 v = a->a + a->i[i] + shift; 791 n = a->i[i+1] - a->i[i]; 792 alpha = x[i]; 793 while (n-->0) {y[*idx++] += alpha * *v++;} 794 } 795 PetscLogFlops(2*a->nz - A->n); 796 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 797 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 798 PetscFunctionReturn(0); 799 } 800 801 #undef __FUNCT__ 802 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 803 int MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 804 { 805 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 806 Scalar *x,*y,*v,alpha; 807 int ierr,m = A->m,n,i,*idx,shift = a->indexshift; 808 809 PetscFunctionBegin; 810 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 811 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 812 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 813 y = y + shift; /* shift for Fortran start by 1 indexing */ 814 for (i=0; i<m; i++) { 815 idx = a->j + a->i[i] + shift; 816 v = a->a + a->i[i] + shift; 817 n = a->i[i+1] - a->i[i]; 818 alpha = x[i]; 819 while (n-->0) {y[*idx++] += alpha * *v++;} 820 } 821 PetscLogFlops(2*a->nz); 822 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 823 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 824 PetscFunctionReturn(0); 825 } 826 827 #undef __FUNCT__ 828 #define __FUNCT__ "MatMult_SeqAIJ" 829 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 830 { 831 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 832 Scalar *x,*y,*v,sum; 833 int ierr,m = A->m,*idx,shift = a->indexshift,*ii; 834 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 835 int n,i,jrow,j; 836 #endif 837 838 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 839 #pragma disjoint(*x,*y,*v) 840 #endif 841 842 PetscFunctionBegin; 843 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 844 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 845 x = x + shift; /* shift for Fortran start by 1 indexing */ 846 idx = a->j; 847 v = a->a; 848 ii = a->i; 849 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 850 fortranmultaij_(&m,x,ii,idx+shift,v+shift,y); 851 #else 852 v += shift; /* shift for Fortran start by 1 indexing */ 853 idx += shift; 854 for (i=0; i<m; i++) { 855 jrow = ii[i]; 856 n = ii[i+1] - jrow; 857 sum = 0.0; 858 for (j=0; j<n; j++) { 859 sum += v[jrow]*x[idx[jrow]]; jrow++; 860 } 861 y[i] = sum; 862 } 863 #endif 864 PetscLogFlops(2*a->nz - m); 865 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 866 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 867 PetscFunctionReturn(0); 868 } 869 870 #undef __FUNCT__ 871 #define __FUNCT__ "MatMultAdd_SeqAIJ" 872 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 873 { 874 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 875 Scalar *x,*y,*z,*v,sum; 876 int ierr,m = A->m,*idx,shift = a->indexshift,*ii; 877 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 878 int n,i,jrow,j; 879 #endif 880 881 PetscFunctionBegin; 882 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 883 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 884 if (zz != yy) { 885 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 886 } else { 887 z = y; 888 } 889 x = x + shift; /* shift for Fortran start by 1 indexing */ 890 idx = a->j; 891 v = a->a; 892 ii = a->i; 893 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 894 fortranmultaddaij_(&m,x,ii,idx+shift,v+shift,y,z); 895 #else 896 v += shift; /* shift for Fortran start by 1 indexing */ 897 idx += shift; 898 for (i=0; i<m; i++) { 899 jrow = ii[i]; 900 n = ii[i+1] - jrow; 901 sum = y[i]; 902 for (j=0; j<n; j++) { 903 sum += v[jrow]*x[idx[jrow]]; jrow++; 904 } 905 z[i] = sum; 906 } 907 #endif 908 PetscLogFlops(2*a->nz); 909 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 910 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 911 if (zz != yy) { 912 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 913 } 914 PetscFunctionReturn(0); 915 } 916 917 /* 918 Adds diagonal pointers to sparse matrix structure. 919 */ 920 #undef __FUNCT__ 921 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 922 int MatMarkDiagonal_SeqAIJ(Mat A) 923 { 924 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 925 int i,j,*diag,m = A->m,shift = a->indexshift,ierr; 926 927 PetscFunctionBegin; 928 if (a->diag) PetscFunctionReturn(0); 929 930 ierr = PetscMalloc((m+1)*sizeof(int),&diag);CHKERRQ(ierr); 931 PetscLogObjectMemory(A,(m+1)*sizeof(int)); 932 for (i=0; i<A->m; i++) { 933 diag[i] = a->i[i+1]; 934 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 935 if (a->j[j]+shift == i) { 936 diag[i] = j - shift; 937 break; 938 } 939 } 940 } 941 a->diag = diag; 942 PetscFunctionReturn(0); 943 } 944 945 /* 946 Checks for missing diagonals 947 */ 948 #undef __FUNCT__ 949 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 950 int MatMissingDiagonal_SeqAIJ(Mat A) 951 { 952 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 953 int *diag,*jj = a->j,i,shift = a->indexshift,ierr; 954 955 PetscFunctionBegin; 956 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 957 diag = a->diag; 958 for (i=0; i<A->m; i++) { 959 if (jj[diag[i]+shift] != i-shift) { 960 SETERRQ1(1,"Matrix is missing diagonal number %d",i); 961 } 962 } 963 PetscFunctionReturn(0); 964 } 965 966 #undef __FUNCT__ 967 #define __FUNCT__ "MatRelax_SeqAIJ" 968 int MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,int its,Vec xx) 969 { 970 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 971 Scalar *x,*b,*bs, d,*xs,sum,*v = a->a,*t=0,scale,*ts,*xb,*idiag=0; 972 int ierr,*idx,*diag,n = A->n,m = A->m,i,shift = a->indexshift; 973 974 PetscFunctionBegin; 975 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 976 if (xx != bb) { 977 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 978 } else { 979 b = x; 980 } 981 982 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 983 diag = a->diag; 984 xs = x + shift; /* shifted by one for index start of a or a->j*/ 985 if (flag == SOR_APPLY_UPPER) { 986 /* apply (U + D/omega) to the vector */ 987 bs = b + shift; 988 for (i=0; i<m; i++) { 989 d = fshift + a->a[diag[i] + shift]; 990 n = a->i[i+1] - diag[i] - 1; 991 PetscLogFlops(2*n-1); 992 idx = a->j + diag[i] + (!shift); 993 v = a->a + diag[i] + (!shift); 994 sum = b[i]*d/omega; 995 SPARSEDENSEDOT(sum,bs,v,idx,n); 996 x[i] = sum; 997 } 998 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 999 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1000 PetscFunctionReturn(0); 1001 } 1002 1003 /* setup workspace for Eisenstat */ 1004 if (flag & SOR_EISENSTAT) { 1005 if (!a->idiag) { 1006 ierr = PetscMalloc(2*m*sizeof(Scalar),&a->idiag);CHKERRQ(ierr); 1007 a->ssor = a->idiag + m; 1008 v = a->a; 1009 for (i=0; i<m; i++) { a->idiag[i] = 1.0/v[diag[i]];} 1010 } 1011 t = a->ssor; 1012 idiag = a->idiag; 1013 } 1014 /* Let A = L + U + D; where L is lower trianglar, 1015 U is upper triangular, E is diagonal; This routine applies 1016 1017 (L + E)^{-1} A (U + E)^{-1} 1018 1019 to a vector efficiently using Eisenstat's trick. This is for 1020 the case of SSOR preconditioner, so E is D/omega where omega 1021 is the relaxation factor. 1022 */ 1023 1024 if (flag == SOR_APPLY_LOWER) { 1025 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1026 } else if ((flag & SOR_EISENSTAT) && omega == 1.0 && shift == 0 && fshift == 0.0) { 1027 /* special case for omega = 1.0 saves flops and some integer ops */ 1028 Scalar *v2; 1029 1030 v2 = a->a; 1031 /* x = (E + U)^{-1} b */ 1032 for (i=m-1; i>=0; i--) { 1033 n = a->i[i+1] - diag[i] - 1; 1034 idx = a->j + diag[i] + 1; 1035 v = a->a + diag[i] + 1; 1036 sum = b[i]; 1037 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1038 x[i] = sum*idiag[i]; 1039 1040 /* t = b - (2*E - D)x */ 1041 t[i] = b[i] - (v2[diag[i]])*x[i]; 1042 } 1043 1044 /* t = (E + L)^{-1}t */ 1045 diag = a->diag; 1046 for (i=0; i<m; i++) { 1047 n = diag[i] - a->i[i]; 1048 idx = a->j + a->i[i]; 1049 v = a->a + a->i[i]; 1050 sum = t[i]; 1051 SPARSEDENSEMDOT(sum,t,v,idx,n); 1052 t[i] = sum*idiag[i]; 1053 1054 /* x = x + t */ 1055 x[i] += t[i]; 1056 } 1057 1058 PetscLogFlops(3*m-1 + 2*a->nz); 1059 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1060 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1061 PetscFunctionReturn(0); 1062 } else if (flag & SOR_EISENSTAT) { 1063 /* Let A = L + U + D; where L is lower trianglar, 1064 U is upper triangular, E is diagonal; This routine applies 1065 1066 (L + E)^{-1} A (U + E)^{-1} 1067 1068 to a vector efficiently using Eisenstat's trick. This is for 1069 the case of SSOR preconditioner, so E is D/omega where omega 1070 is the relaxation factor. 1071 */ 1072 scale = (2.0/omega) - 1.0; 1073 1074 /* x = (E + U)^{-1} b */ 1075 for (i=m-1; i>=0; i--) { 1076 d = fshift + a->a[diag[i] + shift]; 1077 n = a->i[i+1] - diag[i] - 1; 1078 idx = a->j + diag[i] + (!shift); 1079 v = a->a + diag[i] + (!shift); 1080 sum = b[i]; 1081 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1082 x[i] = omega*(sum/d); 1083 } 1084 1085 /* t = b - (2*E - D)x */ 1086 v = a->a; 1087 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 1088 1089 /* t = (E + L)^{-1}t */ 1090 ts = t + shift; /* shifted by one for index start of a or a->j*/ 1091 diag = a->diag; 1092 for (i=0; i<m; i++) { 1093 d = fshift + a->a[diag[i]+shift]; 1094 n = diag[i] - a->i[i]; 1095 idx = a->j + a->i[i] + shift; 1096 v = a->a + a->i[i] + shift; 1097 sum = t[i]; 1098 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1099 t[i] = omega*(sum/d); 1100 /* x = x + t */ 1101 x[i] += t[i]; 1102 } 1103 1104 PetscLogFlops(6*m-1 + 2*a->nz); 1105 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1106 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1107 PetscFunctionReturn(0); 1108 } 1109 if (flag & SOR_ZERO_INITIAL_GUESS) { 1110 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1111 for (i=0; i<m; i++) { 1112 d = fshift + a->a[diag[i]+shift]; 1113 n = diag[i] - a->i[i]; 1114 PetscLogFlops(2*n-1); 1115 idx = a->j + a->i[i] + shift; 1116 v = a->a + a->i[i] + shift; 1117 sum = b[i]; 1118 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1119 x[i] = omega*(sum/d); 1120 } 1121 xb = x; 1122 } else xb = b; 1123 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1124 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1125 for (i=0; i<m; i++) { 1126 x[i] *= a->a[diag[i]+shift]; 1127 } 1128 PetscLogFlops(m); 1129 } 1130 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1131 for (i=m-1; i>=0; i--) { 1132 d = fshift + a->a[diag[i] + shift]; 1133 n = a->i[i+1] - diag[i] - 1; 1134 PetscLogFlops(2*n-1); 1135 idx = a->j + diag[i] + (!shift); 1136 v = a->a + diag[i] + (!shift); 1137 sum = xb[i]; 1138 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1139 x[i] = omega*(sum/d); 1140 } 1141 } 1142 its--; 1143 } 1144 while (its--) { 1145 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1146 for (i=0; i<m; i++) { 1147 d = fshift + a->a[diag[i]+shift]; 1148 n = a->i[i+1] - a->i[i]; 1149 PetscLogFlops(2*n-1); 1150 idx = a->j + a->i[i] + shift; 1151 v = a->a + a->i[i] + shift; 1152 sum = b[i]; 1153 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1154 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1155 } 1156 } 1157 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1158 for (i=m-1; i>=0; i--) { 1159 d = fshift + a->a[diag[i] + shift]; 1160 n = a->i[i+1] - a->i[i]; 1161 PetscLogFlops(2*n-1); 1162 idx = a->j + a->i[i] + shift; 1163 v = a->a + a->i[i] + shift; 1164 sum = b[i]; 1165 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1166 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1167 } 1168 } 1169 } 1170 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1171 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1172 PetscFunctionReturn(0); 1173 } 1174 1175 #undef __FUNCT__ 1176 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1177 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1178 { 1179 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1180 1181 PetscFunctionBegin; 1182 info->rows_global = (double)A->m; 1183 info->columns_global = (double)A->n; 1184 info->rows_local = (double)A->m; 1185 info->columns_local = (double)A->n; 1186 info->block_size = 1.0; 1187 info->nz_allocated = (double)a->maxnz; 1188 info->nz_used = (double)a->nz; 1189 info->nz_unneeded = (double)(a->maxnz - a->nz); 1190 info->assemblies = (double)A->num_ass; 1191 info->mallocs = (double)a->reallocs; 1192 info->memory = A->mem; 1193 if (A->factor) { 1194 info->fill_ratio_given = A->info.fill_ratio_given; 1195 info->fill_ratio_needed = A->info.fill_ratio_needed; 1196 info->factor_mallocs = A->info.factor_mallocs; 1197 } else { 1198 info->fill_ratio_given = 0; 1199 info->fill_ratio_needed = 0; 1200 info->factor_mallocs = 0; 1201 } 1202 PetscFunctionReturn(0); 1203 } 1204 1205 EXTERN int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,MatLUInfo*,Mat*); 1206 EXTERN int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1207 EXTERN int MatLUFactor_SeqAIJ(Mat,IS,IS,MatLUInfo*); 1208 EXTERN int MatSolve_SeqAIJ(Mat,Vec,Vec); 1209 EXTERN int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1210 EXTERN int MatSolveTranspose_SeqAIJ(Mat,Vec,Vec); 1211 EXTERN int MatSolveTransposeAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1212 1213 #undef __FUNCT__ 1214 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1215 int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 1216 { 1217 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1218 int i,ierr,N,*rows,m = A->m - 1,shift = a->indexshift; 1219 1220 PetscFunctionBegin; 1221 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1222 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1223 if (a->keepzeroedrows) { 1224 for (i=0; i<N; i++) { 1225 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1226 ierr = PetscMemzero(&a->a[a->i[rows[i]]+shift],a->ilen[rows[i]]*sizeof(Scalar));CHKERRQ(ierr); 1227 } 1228 if (diag) { 1229 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1230 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1231 for (i=0; i<N; i++) { 1232 a->a[a->diag[rows[i]]] = *diag; 1233 } 1234 } 1235 } else { 1236 if (diag) { 1237 for (i=0; i<N; i++) { 1238 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1239 if (a->ilen[rows[i]] > 0) { 1240 a->ilen[rows[i]] = 1; 1241 a->a[a->i[rows[i]]+shift] = *diag; 1242 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1243 } else { /* in case row was completely empty */ 1244 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1245 } 1246 } 1247 } else { 1248 for (i=0; i<N; i++) { 1249 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1250 a->ilen[rows[i]] = 0; 1251 } 1252 } 1253 } 1254 ierr = ISRestoreIndices(is,&rows);CHKERRQ(ierr); 1255 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1256 PetscFunctionReturn(0); 1257 } 1258 1259 #undef __FUNCT__ 1260 #define __FUNCT__ "MatGetOwnershipRange_SeqAIJ" 1261 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1262 { 1263 PetscFunctionBegin; 1264 if (m) *m = 0; 1265 if (n) *n = A->m; 1266 PetscFunctionReturn(0); 1267 } 1268 1269 #undef __FUNCT__ 1270 #define __FUNCT__ "MatGetRow_SeqAIJ" 1271 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1272 { 1273 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1274 int *itmp,i,shift = a->indexshift,ierr; 1275 1276 PetscFunctionBegin; 1277 if (row < 0 || row >= A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %d out of range",row); 1278 1279 *nz = a->i[row+1] - a->i[row]; 1280 if (v) *v = a->a + a->i[row] + shift; 1281 if (idx) { 1282 itmp = a->j + a->i[row] + shift; 1283 if (*nz && shift) { 1284 ierr = PetscMalloc((*nz)*sizeof(int),idx);CHKERRQ(ierr); 1285 for (i=0; i<(*nz); i++) {(*idx)[i] = itmp[i] + shift;} 1286 } else if (*nz) { 1287 *idx = itmp; 1288 } 1289 else *idx = 0; 1290 } 1291 PetscFunctionReturn(0); 1292 } 1293 1294 #undef __FUNCT__ 1295 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1296 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1297 { 1298 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1299 int ierr; 1300 1301 PetscFunctionBegin; 1302 if (idx) {if (*idx && a->indexshift) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} 1303 PetscFunctionReturn(0); 1304 } 1305 1306 #undef __FUNCT__ 1307 #define __FUNCT__ "MatNorm_SeqAIJ" 1308 int MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *norm) 1309 { 1310 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1311 Scalar *v = a->a; 1312 PetscReal sum = 0.0; 1313 int i,j,shift = a->indexshift,ierr; 1314 1315 PetscFunctionBegin; 1316 if (type == NORM_FROBENIUS) { 1317 for (i=0; i<a->nz; i++) { 1318 #if defined(PETSC_USE_COMPLEX) 1319 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1320 #else 1321 sum += (*v)*(*v); v++; 1322 #endif 1323 } 1324 *norm = sqrt(sum); 1325 } else if (type == NORM_1) { 1326 PetscReal *tmp; 1327 int *jj = a->j; 1328 ierr = PetscMalloc((A->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1329 ierr = PetscMemzero(tmp,A->n*sizeof(PetscReal));CHKERRQ(ierr); 1330 *norm = 0.0; 1331 for (j=0; j<a->nz; j++) { 1332 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1333 } 1334 for (j=0; j<A->n; j++) { 1335 if (tmp[j] > *norm) *norm = tmp[j]; 1336 } 1337 ierr = PetscFree(tmp);CHKERRQ(ierr); 1338 } else if (type == NORM_INFINITY) { 1339 *norm = 0.0; 1340 for (j=0; j<A->m; j++) { 1341 v = a->a + a->i[j] + shift; 1342 sum = 0.0; 1343 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1344 sum += PetscAbsScalar(*v); v++; 1345 } 1346 if (sum > *norm) *norm = sum; 1347 } 1348 } else { 1349 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1350 } 1351 PetscFunctionReturn(0); 1352 } 1353 1354 #undef __FUNCT__ 1355 #define __FUNCT__ "MatTranspose_SeqAIJ" 1356 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1357 { 1358 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1359 Mat C; 1360 int i,ierr,*aj = a->j,*ai = a->i,m = A->m,len,*col; 1361 int shift = a->indexshift; 1362 Scalar *array = a->a; 1363 1364 PetscFunctionBegin; 1365 if (!B && m != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1366 ierr = PetscMalloc((1+A->n)*sizeof(int),&col);CHKERRQ(ierr); 1367 ierr = PetscMemzero(col,(1+A->n)*sizeof(int));CHKERRQ(ierr); 1368 if (shift) { 1369 for (i=0; i<ai[m]-1; i++) aj[i] -= 1; 1370 } 1371 for (i=0; i<ai[m]+shift; i++) col[aj[i]] += 1; 1372 ierr = MatCreateSeqAIJ(A->comm,A->n,m,0,col,&C);CHKERRQ(ierr); 1373 ierr = PetscFree(col);CHKERRQ(ierr); 1374 for (i=0; i<m; i++) { 1375 len = ai[i+1]-ai[i]; 1376 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1377 array += len; 1378 aj += len; 1379 } 1380 if (shift) { 1381 for (i=0; i<ai[m]-1; i++) aj[i] += 1; 1382 } 1383 1384 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1385 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1386 1387 if (B) { 1388 *B = C; 1389 } else { 1390 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1391 } 1392 PetscFunctionReturn(0); 1393 } 1394 1395 #undef __FUNCT__ 1396 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1397 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1398 { 1399 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1400 Scalar *l,*r,x,*v; 1401 int ierr,i,j,m = A->m,n = A->n,M,nz = a->nz,*jj,shift = a->indexshift; 1402 1403 PetscFunctionBegin; 1404 if (ll) { 1405 /* The local size is used so that VecMPI can be passed to this routine 1406 by MatDiagonalScale_MPIAIJ */ 1407 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1408 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1409 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1410 v = a->a; 1411 for (i=0; i<m; i++) { 1412 x = l[i]; 1413 M = a->i[i+1] - a->i[i]; 1414 for (j=0; j<M; j++) { (*v++) *= x;} 1415 } 1416 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1417 PetscLogFlops(nz); 1418 } 1419 if (rr) { 1420 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1421 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1422 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1423 v = a->a; jj = a->j; 1424 for (i=0; i<nz; i++) { 1425 (*v++) *= r[*jj++ + shift]; 1426 } 1427 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1428 PetscLogFlops(nz); 1429 } 1430 PetscFunctionReturn(0); 1431 } 1432 1433 #undef __FUNCT__ 1434 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1435 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatReuse scall,Mat *B) 1436 { 1437 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1438 int *smap,i,k,kstart,kend,ierr,oldcols = A->n,*lens; 1439 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1440 int *irow,*icol,nrows,ncols,shift = a->indexshift,*ssmap; 1441 int *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1442 Scalar *a_new,*mat_a; 1443 Mat C; 1444 PetscTruth stride; 1445 1446 PetscFunctionBegin; 1447 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1448 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1449 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1450 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1451 1452 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1453 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1454 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1455 1456 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1457 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1458 if (stride && step == 1) { 1459 /* special case of contiguous rows */ 1460 ierr = PetscMalloc((ncols+nrows+1)*sizeof(int),&lens);CHKERRQ(ierr); 1461 starts = lens + ncols; 1462 /* loop over new rows determining lens and starting points */ 1463 for (i=0; i<nrows; i++) { 1464 kstart = ai[irow[i]]+shift; 1465 kend = kstart + ailen[irow[i]]; 1466 for (k=kstart; k<kend; k++) { 1467 if (aj[k]+shift >= first) { 1468 starts[i] = k; 1469 break; 1470 } 1471 } 1472 sum = 0; 1473 while (k < kend) { 1474 if (aj[k++]+shift >= first+ncols) break; 1475 sum++; 1476 } 1477 lens[i] = sum; 1478 } 1479 /* create submatrix */ 1480 if (scall == MAT_REUSE_MATRIX) { 1481 int n_cols,n_rows; 1482 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1483 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1484 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1485 C = *B; 1486 } else { 1487 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1488 } 1489 c = (Mat_SeqAIJ*)C->data; 1490 1491 /* loop over rows inserting into submatrix */ 1492 a_new = c->a; 1493 j_new = c->j; 1494 i_new = c->i; 1495 i_new[0] = -shift; 1496 for (i=0; i<nrows; i++) { 1497 ii = starts[i]; 1498 lensi = lens[i]; 1499 for (k=0; k<lensi; k++) { 1500 *j_new++ = aj[ii+k] - first; 1501 } 1502 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar));CHKERRQ(ierr); 1503 a_new += lensi; 1504 i_new[i+1] = i_new[i] + lensi; 1505 c->ilen[i] = lensi; 1506 } 1507 ierr = PetscFree(lens);CHKERRQ(ierr); 1508 } else { 1509 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1510 ierr = PetscMalloc((1+oldcols)*sizeof(int),&smap);CHKERRQ(ierr); 1511 ssmap = smap + shift; 1512 ierr = PetscMalloc((1+nrows)*sizeof(int),&lens);CHKERRQ(ierr); 1513 ierr = PetscMemzero(smap,oldcols*sizeof(int));CHKERRQ(ierr); 1514 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1515 /* determine lens of each row */ 1516 for (i=0; i<nrows; i++) { 1517 kstart = ai[irow[i]]+shift; 1518 kend = kstart + a->ilen[irow[i]]; 1519 lens[i] = 0; 1520 for (k=kstart; k<kend; k++) { 1521 if (ssmap[aj[k]]) { 1522 lens[i]++; 1523 } 1524 } 1525 } 1526 /* Create and fill new matrix */ 1527 if (scall == MAT_REUSE_MATRIX) { 1528 PetscTruth equal; 1529 1530 c = (Mat_SeqAIJ *)((*B)->data); 1531 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1532 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(int),&equal);CHKERRQ(ierr); 1533 if (!equal) { 1534 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1535 } 1536 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(int));CHKERRQ(ierr); 1537 C = *B; 1538 } else { 1539 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1540 } 1541 c = (Mat_SeqAIJ *)(C->data); 1542 for (i=0; i<nrows; i++) { 1543 row = irow[i]; 1544 kstart = ai[row]+shift; 1545 kend = kstart + a->ilen[row]; 1546 mat_i = c->i[i]+shift; 1547 mat_j = c->j + mat_i; 1548 mat_a = c->a + mat_i; 1549 mat_ilen = c->ilen + i; 1550 for (k=kstart; k<kend; k++) { 1551 if ((tcol=ssmap[a->j[k]])) { 1552 *mat_j++ = tcol - (!shift); 1553 *mat_a++ = a->a[k]; 1554 (*mat_ilen)++; 1555 1556 } 1557 } 1558 } 1559 /* Free work space */ 1560 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1561 ierr = PetscFree(smap);CHKERRQ(ierr); 1562 ierr = PetscFree(lens);CHKERRQ(ierr); 1563 } 1564 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1565 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1566 1567 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1568 *B = C; 1569 PetscFunctionReturn(0); 1570 } 1571 1572 /* 1573 */ 1574 #undef __FUNCT__ 1575 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1576 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatILUInfo *info) 1577 { 1578 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1579 int ierr; 1580 Mat outA; 1581 PetscTruth row_identity,col_identity; 1582 1583 PetscFunctionBegin; 1584 if (info && info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1585 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1586 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1587 if (!row_identity || !col_identity) { 1588 SETERRQ(1,"Row and column permutations must be identity for in-place ILU"); 1589 } 1590 1591 outA = inA; 1592 inA->factor = FACTOR_LU; 1593 a->row = row; 1594 a->col = col; 1595 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1596 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1597 1598 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1599 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1600 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1601 PetscLogObjectParent(inA,a->icol); 1602 1603 if (!a->solve_work) { /* this matrix may have been factored before */ 1604 ierr = PetscMalloc((inA->m+1)*sizeof(Scalar),&a->solve_work);CHKERRQ(ierr); 1605 } 1606 1607 if (!a->diag) { 1608 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1609 } 1610 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1611 PetscFunctionReturn(0); 1612 } 1613 1614 #include "petscblaslapack.h" 1615 #undef __FUNCT__ 1616 #define __FUNCT__ "MatScale_SeqAIJ" 1617 int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1618 { 1619 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1620 int one = 1; 1621 1622 PetscFunctionBegin; 1623 BLscal_(&a->nz,alpha,a->a,&one); 1624 PetscLogFlops(a->nz); 1625 PetscFunctionReturn(0); 1626 } 1627 1628 #undef __FUNCT__ 1629 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1630 int MatGetSubMatrices_SeqAIJ(Mat A,int n,IS *irow,IS *icol,MatReuse scall,Mat **B) 1631 { 1632 int ierr,i; 1633 1634 PetscFunctionBegin; 1635 if (scall == MAT_INITIAL_MATRIX) { 1636 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1637 } 1638 1639 for (i=0; i<n; i++) { 1640 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1641 } 1642 PetscFunctionReturn(0); 1643 } 1644 1645 #undef __FUNCT__ 1646 #define __FUNCT__ "MatGetBlockSize_SeqAIJ" 1647 int MatGetBlockSize_SeqAIJ(Mat A,int *bs) 1648 { 1649 PetscFunctionBegin; 1650 *bs = 1; 1651 PetscFunctionReturn(0); 1652 } 1653 1654 #undef __FUNCT__ 1655 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1656 int MatIncreaseOverlap_SeqAIJ(Mat A,int is_max,IS *is,int ov) 1657 { 1658 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1659 int shift,row,i,j,k,l,m,n,*idx,ierr,*nidx,isz,val; 1660 int start,end,*ai,*aj; 1661 PetscBT table; 1662 1663 PetscFunctionBegin; 1664 shift = a->indexshift; 1665 m = A->m; 1666 ai = a->i; 1667 aj = a->j+shift; 1668 1669 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal overlap value used"); 1670 1671 ierr = PetscMalloc((m+1)*sizeof(int),&nidx);CHKERRQ(ierr); 1672 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1673 1674 for (i=0; i<is_max; i++) { 1675 /* Initialize the two local arrays */ 1676 isz = 0; 1677 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1678 1679 /* Extract the indices, assume there can be duplicate entries */ 1680 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1681 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1682 1683 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1684 for (j=0; j<n ; ++j){ 1685 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1686 } 1687 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1688 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1689 1690 k = 0; 1691 for (j=0; j<ov; j++){ /* for each overlap */ 1692 n = isz; 1693 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1694 row = nidx[k]; 1695 start = ai[row]; 1696 end = ai[row+1]; 1697 for (l = start; l<end ; l++){ 1698 val = aj[l] + shift; 1699 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1700 } 1701 } 1702 } 1703 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1704 } 1705 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1706 ierr = PetscFree(nidx);CHKERRQ(ierr); 1707 PetscFunctionReturn(0); 1708 } 1709 1710 /* -------------------------------------------------------------- */ 1711 #undef __FUNCT__ 1712 #define __FUNCT__ "MatPermute_SeqAIJ" 1713 int MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1714 { 1715 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1716 Scalar *vwork; 1717 int i,ierr,nz,m = A->m,n = A->n,*cwork; 1718 int *row,*col,*cnew,j,*lens; 1719 IS icolp,irowp; 1720 1721 PetscFunctionBegin; 1722 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1723 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1724 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1725 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1726 1727 /* determine lengths of permuted rows */ 1728 ierr = PetscMalloc((m+1)*sizeof(int),&lens);CHKERRQ(ierr); 1729 for (i=0; i<m; i++) { 1730 lens[row[i]] = a->i[i+1] - a->i[i]; 1731 } 1732 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1733 ierr = PetscFree(lens);CHKERRQ(ierr); 1734 1735 ierr = PetscMalloc(n*sizeof(int),&cnew);CHKERRQ(ierr); 1736 for (i=0; i<m; i++) { 1737 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1738 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1739 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1740 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1741 } 1742 ierr = PetscFree(cnew);CHKERRQ(ierr); 1743 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1744 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1745 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1746 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1747 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1748 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1749 PetscFunctionReturn(0); 1750 } 1751 1752 #undef __FUNCT__ 1753 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1754 int MatPrintHelp_SeqAIJ(Mat A) 1755 { 1756 static PetscTruth called = PETSC_FALSE; 1757 MPI_Comm comm = A->comm; 1758 int ierr; 1759 1760 PetscFunctionBegin; 1761 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1762 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1763 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1764 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1765 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1766 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1767 #if defined(PETSC_HAVE_ESSL) 1768 ierr = (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n");CHKERRQ(ierr); 1769 #endif 1770 #if defined(PETSC_HAVE_LUSOL) 1771 ierr = (*PetscHelpPrintf)(comm," -mat_aij_lusol: Use the Stanford LUSOL sparse factorization and solve.\n");CHKERRQ(ierr); 1772 #endif 1773 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) 1774 ierr = (*PetscHelpPrintf)(comm," -mat_aij_matlab: Use Matlab engine sparse LU factorization and solve.\n");CHKERRQ(ierr); 1775 #endif 1776 PetscFunctionReturn(0); 1777 } 1778 EXTERN int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg); 1779 EXTERN int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1780 EXTERN int MatILUDTFactor_SeqAIJ(Mat,MatILUInfo*,IS,IS,Mat*); 1781 #undef __FUNCT__ 1782 #define __FUNCT__ "MatCopy_SeqAIJ" 1783 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1784 { 1785 int ierr; 1786 PetscTruth flg; 1787 1788 PetscFunctionBegin; 1789 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flg);CHKERRQ(ierr); 1790 if (str == SAME_NONZERO_PATTERN && flg) { 1791 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1792 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1793 1794 if (a->i[A->m]+a->indexshift != b->i[B->m]+a->indexshift) { 1795 SETERRQ(1,"Number of nonzeros in two matrices are different"); 1796 } 1797 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m]+a->indexshift)*sizeof(Scalar));CHKERRQ(ierr); 1798 } else { 1799 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1800 } 1801 PetscFunctionReturn(0); 1802 } 1803 1804 #undef __FUNCT__ 1805 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1806 int MatSetUpPreallocation_SeqAIJ(Mat A) 1807 { 1808 int ierr; 1809 1810 PetscFunctionBegin; 1811 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1812 PetscFunctionReturn(0); 1813 } 1814 1815 #undef __FUNCT__ 1816 #define __FUNCT__ "MatGetArray_SeqAIJ" 1817 int MatGetArray_SeqAIJ(Mat A,Scalar **array) 1818 { 1819 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1820 PetscFunctionBegin; 1821 *array = a->a; 1822 PetscFunctionReturn(0); 1823 } 1824 1825 #undef __FUNCT__ 1826 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1827 int MatRestoreArray_SeqAIJ(Mat A,Scalar **array) 1828 { 1829 PetscFunctionBegin; 1830 PetscFunctionReturn(0); 1831 } 1832 1833 #undef __FUNCT__ 1834 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1835 int MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1836 { 1837 int (*f)(void *,Vec,Vec,void*) = (int (*)(void *,Vec,Vec,void *))coloring->f; 1838 int k,ierr,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1839 Scalar dx,mone = -1.0,*y,*xx,*w3_array; 1840 Scalar *vscale_array; 1841 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1842 Vec w1,w2,w3; 1843 void *fctx = coloring->fctx; 1844 PetscTruth flg; 1845 1846 PetscFunctionBegin; 1847 if (!coloring->w1) { 1848 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1849 PetscLogObjectParent(coloring,coloring->w1); 1850 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1851 PetscLogObjectParent(coloring,coloring->w2); 1852 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1853 PetscLogObjectParent(coloring,coloring->w3); 1854 } 1855 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1856 1857 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1858 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1859 if (flg) { 1860 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 1861 } else { 1862 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1863 } 1864 1865 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1866 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1867 1868 /* 1869 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1870 coloring->F for the coarser grids from the finest 1871 */ 1872 if (coloring->F) { 1873 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1874 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1875 if (m1 != m2) { 1876 coloring->F = 0; 1877 } 1878 } 1879 1880 if (coloring->F) { 1881 w1 = coloring->F; 1882 coloring->F = 0; 1883 } else { 1884 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1885 } 1886 1887 /* 1888 Compute all the scale factors and share with other processors 1889 */ 1890 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 1891 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 1892 for (k=0; k<coloring->ncolors; k++) { 1893 /* 1894 Loop over each column associated with color adding the 1895 perturbation to the vector w3. 1896 */ 1897 for (l=0; l<coloring->ncolumns[k]; l++) { 1898 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1899 dx = xx[col]; 1900 if (dx == 0.0) dx = 1.0; 1901 #if !defined(PETSC_USE_COMPLEX) 1902 if (dx < umin && dx >= 0.0) dx = umin; 1903 else if (dx < 0.0 && dx > -umin) dx = -umin; 1904 #else 1905 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1906 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1907 #endif 1908 dx *= epsilon; 1909 vscale_array[col] = 1.0/dx; 1910 } 1911 } 1912 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1913 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1914 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1915 1916 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 1917 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 1918 1919 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 1920 else vscaleforrow = coloring->columnsforrow; 1921 1922 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1923 /* 1924 Loop over each color 1925 */ 1926 for (k=0; k<coloring->ncolors; k++) { 1927 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 1928 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 1929 /* 1930 Loop over each column associated with color adding the 1931 perturbation to the vector w3. 1932 */ 1933 for (l=0; l<coloring->ncolumns[k]; l++) { 1934 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1935 dx = xx[col]; 1936 if (dx == 0.0) dx = 1.0; 1937 #if !defined(PETSC_USE_COMPLEX) 1938 if (dx < umin && dx >= 0.0) dx = umin; 1939 else if (dx < 0.0 && dx > -umin) dx = -umin; 1940 #else 1941 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1942 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1943 #endif 1944 dx *= epsilon; 1945 if (!PetscAbsScalar(dx)) SETERRQ(1,"Computed 0 differencing parameter"); 1946 w3_array[col] += dx; 1947 } 1948 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 1949 1950 /* 1951 Evaluate function at x1 + dx (here dx is a vector of perturbations) 1952 */ 1953 1954 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 1955 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 1956 1957 /* 1958 Loop over rows of vector, putting results into Jacobian matrix 1959 */ 1960 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 1961 for (l=0; l<coloring->nrows[k]; l++) { 1962 row = coloring->rows[k][l]; 1963 col = coloring->columnsforrow[k][l]; 1964 y[row] *= vscale_array[vscaleforrow[k][l]]; 1965 srow = row + start; 1966 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 1967 } 1968 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 1969 } 1970 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1971 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 1972 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1973 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1974 PetscFunctionReturn(0); 1975 } 1976 1977 /* -------------------------------------------------------------------*/ 1978 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 1979 MatGetRow_SeqAIJ, 1980 MatRestoreRow_SeqAIJ, 1981 MatMult_SeqAIJ, 1982 MatMultAdd_SeqAIJ, 1983 MatMultTranspose_SeqAIJ, 1984 MatMultTransposeAdd_SeqAIJ, 1985 MatSolve_SeqAIJ, 1986 MatSolveAdd_SeqAIJ, 1987 MatSolveTranspose_SeqAIJ, 1988 MatSolveTransposeAdd_SeqAIJ, 1989 MatLUFactor_SeqAIJ, 1990 0, 1991 MatRelax_SeqAIJ, 1992 MatTranspose_SeqAIJ, 1993 MatGetInfo_SeqAIJ, 1994 MatEqual_SeqAIJ, 1995 MatGetDiagonal_SeqAIJ, 1996 MatDiagonalScale_SeqAIJ, 1997 MatNorm_SeqAIJ, 1998 0, 1999 MatAssemblyEnd_SeqAIJ, 2000 MatCompress_SeqAIJ, 2001 MatSetOption_SeqAIJ, 2002 MatZeroEntries_SeqAIJ, 2003 MatZeroRows_SeqAIJ, 2004 MatLUFactorSymbolic_SeqAIJ, 2005 MatLUFactorNumeric_SeqAIJ, 2006 0, 2007 0, 2008 MatSetUpPreallocation_SeqAIJ, 2009 0, 2010 MatGetOwnershipRange_SeqAIJ, 2011 MatILUFactorSymbolic_SeqAIJ, 2012 0, 2013 MatGetArray_SeqAIJ, 2014 MatRestoreArray_SeqAIJ, 2015 MatDuplicate_SeqAIJ, 2016 0, 2017 0, 2018 MatILUFactor_SeqAIJ, 2019 0, 2020 0, 2021 MatGetSubMatrices_SeqAIJ, 2022 MatIncreaseOverlap_SeqAIJ, 2023 MatGetValues_SeqAIJ, 2024 MatCopy_SeqAIJ, 2025 MatPrintHelp_SeqAIJ, 2026 MatScale_SeqAIJ, 2027 0, 2028 0, 2029 MatILUDTFactor_SeqAIJ, 2030 MatGetBlockSize_SeqAIJ, 2031 MatGetRowIJ_SeqAIJ, 2032 MatRestoreRowIJ_SeqAIJ, 2033 MatGetColumnIJ_SeqAIJ, 2034 MatRestoreColumnIJ_SeqAIJ, 2035 MatFDColoringCreate_SeqAIJ, 2036 0, 2037 0, 2038 MatPermute_SeqAIJ, 2039 0, 2040 0, 2041 MatDestroy_SeqAIJ, 2042 MatView_SeqAIJ, 2043 MatGetMaps_Petsc, 2044 0, 2045 0, 2046 0, 2047 0, 2048 0, 2049 0, 2050 0, 2051 0, 2052 MatSetColoring_SeqAIJ, 2053 MatSetValuesAdic_SeqAIJ, 2054 MatSetValuesAdifor_SeqAIJ, 2055 MatFDColoringApply_SeqAIJ}; 2056 2057 EXTERN int MatUseSuperLU_SeqAIJ(Mat); 2058 EXTERN int MatUseEssl_SeqAIJ(Mat); 2059 EXTERN int MatUseLUSOL_SeqAIJ(Mat); 2060 EXTERN int MatUseMatlab_SeqAIJ(Mat); 2061 EXTERN int MatUseDXML_SeqAIJ(Mat); 2062 2063 EXTERN_C_BEGIN 2064 #undef __FUNCT__ 2065 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2066 2067 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 2068 { 2069 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2070 int i,nz,n; 2071 2072 PetscFunctionBegin; 2073 2074 nz = aij->maxnz; 2075 n = mat->n; 2076 for (i=0; i<nz; i++) { 2077 aij->j[i] = indices[i]; 2078 } 2079 aij->nz = nz; 2080 for (i=0; i<n; i++) { 2081 aij->ilen[i] = aij->imax[i]; 2082 } 2083 2084 PetscFunctionReturn(0); 2085 } 2086 EXTERN_C_END 2087 2088 #undef __FUNCT__ 2089 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2090 /*@ 2091 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2092 in the matrix. 2093 2094 Input Parameters: 2095 + mat - the SeqAIJ matrix 2096 - indices - the column indices 2097 2098 Level: advanced 2099 2100 Notes: 2101 This can be called if you have precomputed the nonzero structure of the 2102 matrix and want to provide it to the matrix object to improve the performance 2103 of the MatSetValues() operation. 2104 2105 You MUST have set the correct numbers of nonzeros per row in the call to 2106 MatCreateSeqAIJ(). 2107 2108 MUST be called before any calls to MatSetValues(); 2109 2110 The indices should start with zero, not one. 2111 2112 @*/ 2113 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 2114 { 2115 int ierr,(*f)(Mat,int *); 2116 2117 PetscFunctionBegin; 2118 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2119 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)())&f);CHKERRQ(ierr); 2120 if (f) { 2121 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2122 } else { 2123 SETERRQ(1,"Wrong type of matrix to set column indices"); 2124 } 2125 PetscFunctionReturn(0); 2126 } 2127 2128 /* ----------------------------------------------------------------------------------------*/ 2129 2130 EXTERN_C_BEGIN 2131 #undef __FUNCT__ 2132 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2133 int MatStoreValues_SeqAIJ(Mat mat) 2134 { 2135 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2136 int nz = aij->i[mat->m]+aij->indexshift,ierr; 2137 2138 PetscFunctionBegin; 2139 if (aij->nonew != 1) { 2140 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2141 } 2142 2143 /* allocate space for values if not already there */ 2144 if (!aij->saved_values) { 2145 ierr = PetscMalloc((nz+1)*sizeof(Scalar),&aij->saved_values);CHKERRQ(ierr); 2146 } 2147 2148 /* copy values over */ 2149 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar));CHKERRQ(ierr); 2150 PetscFunctionReturn(0); 2151 } 2152 EXTERN_C_END 2153 2154 #undef __FUNCT__ 2155 #define __FUNCT__ "MatStoreValues" 2156 /*@ 2157 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2158 example, reuse of the linear part of a Jacobian, while recomputing the 2159 nonlinear portion. 2160 2161 Collect on Mat 2162 2163 Input Parameters: 2164 . mat - the matrix (currently on AIJ matrices support this option) 2165 2166 Level: advanced 2167 2168 Common Usage, with SNESSolve(): 2169 $ Create Jacobian matrix 2170 $ Set linear terms into matrix 2171 $ Apply boundary conditions to matrix, at this time matrix must have 2172 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2173 $ boundary conditions again will not change the nonzero structure 2174 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2175 $ ierr = MatStoreValues(mat); 2176 $ Call SNESSetJacobian() with matrix 2177 $ In your Jacobian routine 2178 $ ierr = MatRetrieveValues(mat); 2179 $ Set nonlinear terms in matrix 2180 2181 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2182 $ // build linear portion of Jacobian 2183 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2184 $ ierr = MatStoreValues(mat); 2185 $ loop over nonlinear iterations 2186 $ ierr = MatRetrieveValues(mat); 2187 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2188 $ // call MatAssemblyBegin/End() on matrix 2189 $ Solve linear system with Jacobian 2190 $ endloop 2191 2192 Notes: 2193 Matrix must already be assemblied before calling this routine 2194 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2195 calling this routine. 2196 2197 .seealso: MatRetrieveValues() 2198 2199 @*/ 2200 int MatStoreValues(Mat mat) 2201 { 2202 int ierr,(*f)(Mat); 2203 2204 PetscFunctionBegin; 2205 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2206 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2207 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2208 2209 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)())&f);CHKERRQ(ierr); 2210 if (f) { 2211 ierr = (*f)(mat);CHKERRQ(ierr); 2212 } else { 2213 SETERRQ(1,"Wrong type of matrix to store values"); 2214 } 2215 PetscFunctionReturn(0); 2216 } 2217 2218 EXTERN_C_BEGIN 2219 #undef __FUNCT__ 2220 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2221 int MatRetrieveValues_SeqAIJ(Mat mat) 2222 { 2223 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2224 int nz = aij->i[mat->m]+aij->indexshift,ierr; 2225 2226 PetscFunctionBegin; 2227 if (aij->nonew != 1) { 2228 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2229 } 2230 if (!aij->saved_values) { 2231 SETERRQ(1,"Must call MatStoreValues(A);first"); 2232 } 2233 2234 /* copy values over */ 2235 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(Scalar));CHKERRQ(ierr); 2236 PetscFunctionReturn(0); 2237 } 2238 EXTERN_C_END 2239 2240 #undef __FUNCT__ 2241 #define __FUNCT__ "MatRetrieveValues" 2242 /*@ 2243 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2244 example, reuse of the linear part of a Jacobian, while recomputing the 2245 nonlinear portion. 2246 2247 Collect on Mat 2248 2249 Input Parameters: 2250 . mat - the matrix (currently on AIJ matrices support this option) 2251 2252 Level: advanced 2253 2254 .seealso: MatStoreValues() 2255 2256 @*/ 2257 int MatRetrieveValues(Mat mat) 2258 { 2259 int ierr,(*f)(Mat); 2260 2261 PetscFunctionBegin; 2262 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2263 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2264 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2265 2266 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)())&f);CHKERRQ(ierr); 2267 if (f) { 2268 ierr = (*f)(mat);CHKERRQ(ierr); 2269 } else { 2270 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2271 } 2272 PetscFunctionReturn(0); 2273 } 2274 2275 /* 2276 This allows SeqAIJ matrices to be passed to the matlab engine 2277 */ 2278 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) 2279 #include "engine.h" /* Matlab include file */ 2280 #include "mex.h" /* Matlab include file */ 2281 EXTERN_C_BEGIN 2282 #undef __FUNCT__ 2283 #define __FUNCT__ "MatMatlabEnginePut_SeqAIJ" 2284 int MatMatlabEnginePut_SeqAIJ(PetscObject obj,void *engine) 2285 { 2286 int ierr,i,*ai,*aj; 2287 Mat B = (Mat)obj; 2288 Scalar *array; 2289 mxArray *mat; 2290 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2291 2292 PetscFunctionBegin; 2293 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2294 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(Scalar));CHKERRQ(ierr); 2295 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2296 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2297 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2298 2299 /* Matlab indices start at 0 for sparse (what a surprise) */ 2300 if (aij->indexshift) { 2301 for (i=0; i<B->m+1; i++) { 2302 ai[i]--; 2303 } 2304 for (i=0; i<aij->nz; i++) { 2305 aj[i]--; 2306 } 2307 } 2308 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2309 mxSetName(mat,obj->name); 2310 engPutArray((Engine *)engine,mat); 2311 PetscFunctionReturn(0); 2312 } 2313 EXTERN_C_END 2314 2315 EXTERN_C_BEGIN 2316 #undef __FUNCT__ 2317 #define __FUNCT__ "MatMatlabEngineGet_SeqAIJ" 2318 int MatMatlabEngineGet_SeqAIJ(PetscObject obj,void *engine) 2319 { 2320 int ierr,ii; 2321 Mat mat = (Mat)obj; 2322 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2323 mxArray *mmat; 2324 2325 PetscFunctionBegin; 2326 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2327 aij->indexshift = 0; 2328 2329 mmat = engGetArray((Engine *)engine,obj->name); 2330 2331 aij->nz = (mxGetJc(mmat))[mat->m]; 2332 ierr = PetscMalloc(aij->nz*(sizeof(int)+sizeof(Scalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 2333 aij->j = (int*)(aij->a + aij->nz); 2334 aij->i = aij->j + aij->nz; 2335 aij->singlemalloc = PETSC_TRUE; 2336 aij->freedata = PETSC_TRUE; 2337 2338 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(Scalar));CHKERRQ(ierr); 2339 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2340 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 2341 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 2342 2343 for (ii=0; ii<mat->m; ii++) { 2344 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 2345 } 2346 2347 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2348 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2349 2350 PetscFunctionReturn(0); 2351 } 2352 EXTERN_C_END 2353 #endif 2354 2355 /* --------------------------------------------------------------------------------*/ 2356 2357 #undef __FUNCT__ 2358 #define __FUNCT__ "MatCreateSeqAIJ" 2359 /*@C 2360 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2361 (the default parallel PETSc format). For good matrix assembly performance 2362 the user should preallocate the matrix storage by setting the parameter nz 2363 (or the array nnz). By setting these parameters accurately, performance 2364 during matrix assembly can be increased by more than a factor of 50. 2365 2366 Collective on MPI_Comm 2367 2368 Input Parameters: 2369 + comm - MPI communicator, set to PETSC_COMM_SELF 2370 . m - number of rows 2371 . n - number of columns 2372 . nz - number of nonzeros per row (same for all rows) 2373 - nnz - array containing the number of nonzeros in the various rows 2374 (possibly different for each row) or PETSC_NULL 2375 2376 Output Parameter: 2377 . A - the matrix 2378 2379 Notes: 2380 The AIJ format (also called the Yale sparse matrix format or 2381 compressed row storage), is fully compatible with standard Fortran 77 2382 storage. That is, the stored row and column indices can begin at 2383 either one (as in Fortran) or zero. See the users' manual for details. 2384 2385 Specify the preallocated storage with either nz or nnz (not both). 2386 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2387 allocation. For large problems you MUST preallocate memory or you 2388 will get TERRIBLE performance, see the users' manual chapter on matrices. 2389 2390 By default, this format uses inodes (identical nodes) when possible, to 2391 improve numerical efficiency of matrix-vector products and solves. We 2392 search for consecutive rows with the same nonzero structure, thereby 2393 reusing matrix information to achieve increased efficiency. 2394 2395 Options Database Keys: 2396 + -mat_aij_no_inode - Do not use inodes 2397 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2398 - -mat_aij_oneindex - Internally use indexing starting at 1 2399 rather than 0. Note that when calling MatSetValues(), 2400 the user still MUST index entries starting at 0! 2401 2402 Level: intermediate 2403 2404 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2405 2406 @*/ 2407 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz,Mat *A) 2408 { 2409 int ierr; 2410 2411 PetscFunctionBegin; 2412 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2413 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2414 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2415 PetscFunctionReturn(0); 2416 } 2417 2418 #undef __FUNCT__ 2419 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2420 /*@C 2421 MatSeqAIJSetPreallocation - For good matrix assembly performance 2422 the user should preallocate the matrix storage by setting the parameter nz 2423 (or the array nnz). By setting these parameters accurately, performance 2424 during matrix assembly can be increased by more than a factor of 50. 2425 2426 Collective on MPI_Comm 2427 2428 Input Parameters: 2429 + comm - MPI communicator, set to PETSC_COMM_SELF 2430 . m - number of rows 2431 . n - number of columns 2432 . nz - number of nonzeros per row (same for all rows) 2433 - nnz - array containing the number of nonzeros in the various rows 2434 (possibly different for each row) or PETSC_NULL 2435 2436 Output Parameter: 2437 . A - the matrix 2438 2439 Notes: 2440 The AIJ format (also called the Yale sparse matrix format or 2441 compressed row storage), is fully compatible with standard Fortran 77 2442 storage. That is, the stored row and column indices can begin at 2443 either one (as in Fortran) or zero. See the users' manual for details. 2444 2445 Specify the preallocated storage with either nz or nnz (not both). 2446 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2447 allocation. For large problems you MUST preallocate memory or you 2448 will get TERRIBLE performance, see the users' manual chapter on matrices. 2449 2450 By default, this format uses inodes (identical nodes) when possible, to 2451 improve numerical efficiency of matrix-vector products and solves. We 2452 search for consecutive rows with the same nonzero structure, thereby 2453 reusing matrix information to achieve increased efficiency. 2454 2455 Options Database Keys: 2456 + -mat_aij_no_inode - Do not use inodes 2457 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2458 - -mat_aij_oneindex - Internally use indexing starting at 1 2459 rather than 0. Note that when calling MatSetValues(), 2460 the user still MUST index entries starting at 0! 2461 2462 Level: intermediate 2463 2464 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2465 2466 @*/ 2467 int MatSeqAIJSetPreallocation(Mat B,int nz,int *nnz) 2468 { 2469 Mat_SeqAIJ *b; 2470 int i,len,ierr; 2471 PetscTruth flg2; 2472 2473 PetscFunctionBegin; 2474 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flg2);CHKERRQ(ierr); 2475 if (!flg2) PetscFunctionReturn(0); 2476 2477 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2478 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2479 if (nnz) { 2480 for (i=0; i<B->m; i++) { 2481 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2482 } 2483 } 2484 2485 B->preallocated = PETSC_TRUE; 2486 b = (Mat_SeqAIJ*)B->data; 2487 2488 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2489 if (!nnz) { 2490 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2491 else if (nz <= 0) nz = 1; 2492 for (i=0; i<B->m; i++) b->imax[i] = nz; 2493 nz = nz*B->m; 2494 } else { 2495 nz = 0; 2496 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2497 } 2498 2499 /* allocate the matrix space */ 2500 len = nz*(sizeof(int) + sizeof(Scalar)) + (B->m+1)*sizeof(int); 2501 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2502 b->j = (int*)(b->a + nz); 2503 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2504 b->i = b->j + nz; 2505 b->singlemalloc = PETSC_TRUE; 2506 b->freedata = PETSC_TRUE; 2507 2508 b->i[0] = -b->indexshift; 2509 for (i=1; i<B->m+1; i++) { 2510 b->i[i] = b->i[i-1] + b->imax[i-1]; 2511 } 2512 2513 /* b->ilen will count nonzeros in each row so far. */ 2514 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2515 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2516 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2517 2518 b->nz = 0; 2519 b->maxnz = nz; 2520 B->info.nz_unneeded = (double)b->maxnz; 2521 PetscFunctionReturn(0); 2522 } 2523 2524 EXTERN_C_BEGIN 2525 #undef __FUNCT__ 2526 #define __FUNCT__ "MatCreate_SeqAIJ" 2527 int MatCreate_SeqAIJ(Mat B) 2528 { 2529 Mat_SeqAIJ *b; 2530 int ierr,size; 2531 PetscTruth flg; 2532 2533 PetscFunctionBegin; 2534 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2535 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2536 2537 B->m = B->M = PetscMax(B->m,B->M); 2538 B->n = B->N = PetscMax(B->n,B->N); 2539 2540 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2541 B->data = (void*)b; 2542 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2543 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2544 B->factor = 0; 2545 B->lupivotthreshold = 1.0; 2546 B->mapping = 0; 2547 ierr = PetscOptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2548 ierr = PetscOptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2549 b->row = 0; 2550 b->col = 0; 2551 b->icol = 0; 2552 b->indexshift = 0; 2553 b->reallocs = 0; 2554 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_oneindex",&flg);CHKERRQ(ierr); 2555 if (flg) b->indexshift = -1; 2556 2557 ierr = MapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2558 ierr = MapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2559 2560 b->sorted = PETSC_FALSE; 2561 b->ignorezeroentries = PETSC_FALSE; 2562 b->roworiented = PETSC_TRUE; 2563 b->nonew = 0; 2564 b->diag = 0; 2565 b->solve_work = 0; 2566 b->spptr = 0; 2567 b->inode.use = PETSC_TRUE; 2568 b->inode.node_count = 0; 2569 b->inode.size = 0; 2570 b->inode.limit = 5; 2571 b->inode.max_limit = 5; 2572 b->saved_values = 0; 2573 b->idiag = 0; 2574 b->ssor = 0; 2575 b->keepzeroedrows = PETSC_FALSE; 2576 2577 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2578 2579 /* SuperLU is not currently supported through PETSc */ 2580 #if defined(PETSC_HAVE_SUPERLU) 2581 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_superlu",&flg);CHKERRQ(ierr); 2582 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B);CHKERRQ(ierr); } 2583 #endif 2584 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_essl",&flg);CHKERRQ(ierr); 2585 if (flg) { ierr = MatUseEssl_SeqAIJ(B);CHKERRQ(ierr); } 2586 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_lusol",&flg);CHKERRQ(ierr); 2587 if (flg) { ierr = MatUseLUSOL_SeqAIJ(B);CHKERRQ(ierr); } 2588 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_matlab",&flg);CHKERRQ(ierr); 2589 if (flg) {ierr = MatUseMatlab_SeqAIJ(B);CHKERRQ(ierr);} 2590 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_dxml",&flg);CHKERRQ(ierr); 2591 if (flg) { 2592 if (!b->indexshift) SETERRQ(PETSC_ERR_LIB,"need -mat_aij_oneindex with -mat_aij_dxml"); 2593 ierr = MatUseDXML_SeqAIJ(B);CHKERRQ(ierr); 2594 } 2595 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2596 "MatSeqAIJSetColumnIndices_SeqAIJ", 2597 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2598 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2599 "MatStoreValues_SeqAIJ", 2600 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2601 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2602 "MatRetrieveValues_SeqAIJ", 2603 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2604 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) 2605 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2606 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2607 #endif 2608 PetscFunctionReturn(0); 2609 } 2610 EXTERN_C_END 2611 2612 #undef __FUNCT__ 2613 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2614 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2615 { 2616 Mat C; 2617 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2618 int i,len,m = A->m,shift = a->indexshift,ierr; 2619 2620 PetscFunctionBegin; 2621 *B = 0; 2622 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2623 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 2624 c = (Mat_SeqAIJ*)C->data; 2625 2626 C->factor = A->factor; 2627 c->row = 0; 2628 c->col = 0; 2629 c->icol = 0; 2630 c->indexshift = shift; 2631 c->keepzeroedrows = a->keepzeroedrows; 2632 C->assembled = PETSC_TRUE; 2633 2634 C->M = A->m; 2635 C->N = A->n; 2636 2637 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2638 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2639 for (i=0; i<m; i++) { 2640 c->imax[i] = a->imax[i]; 2641 c->ilen[i] = a->ilen[i]; 2642 } 2643 2644 /* allocate the matrix space */ 2645 c->singlemalloc = PETSC_TRUE; 2646 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 2647 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2648 c->j = (int*)(c->a + a->i[m] + shift); 2649 c->i = c->j + a->i[m] + shift; 2650 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2651 if (m > 0) { 2652 ierr = PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int));CHKERRQ(ierr); 2653 if (cpvalues == MAT_COPY_VALUES) { 2654 ierr = PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar));CHKERRQ(ierr); 2655 } else { 2656 ierr = PetscMemzero(c->a,(a->i[m]+shift)*sizeof(Scalar));CHKERRQ(ierr); 2657 } 2658 } 2659 2660 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2661 c->sorted = a->sorted; 2662 c->roworiented = a->roworiented; 2663 c->nonew = a->nonew; 2664 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2665 c->saved_values = 0; 2666 c->idiag = 0; 2667 c->ssor = 0; 2668 c->ignorezeroentries = a->ignorezeroentries; 2669 c->freedata = PETSC_TRUE; 2670 2671 if (a->diag) { 2672 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2673 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2674 for (i=0; i<m; i++) { 2675 c->diag[i] = a->diag[i]; 2676 } 2677 } else c->diag = 0; 2678 c->inode.use = a->inode.use; 2679 c->inode.limit = a->inode.limit; 2680 c->inode.max_limit = a->inode.max_limit; 2681 if (a->inode.size){ 2682 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2683 c->inode.node_count = a->inode.node_count; 2684 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2685 } else { 2686 c->inode.size = 0; 2687 c->inode.node_count = 0; 2688 } 2689 c->nz = a->nz; 2690 c->maxnz = a->maxnz; 2691 c->solve_work = 0; 2692 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2693 C->preallocated = PETSC_TRUE; 2694 2695 *B = C; 2696 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2697 PetscFunctionReturn(0); 2698 } 2699 2700 EXTERN_C_BEGIN 2701 #undef __FUNCT__ 2702 #define __FUNCT__ "MatLoad_SeqAIJ" 2703 int MatLoad_SeqAIJ(PetscViewer viewer,MatType type,Mat *A) 2704 { 2705 Mat_SeqAIJ *a; 2706 Mat B; 2707 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N,shift; 2708 MPI_Comm comm; 2709 2710 PetscFunctionBegin; 2711 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2712 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2713 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2714 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2715 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2716 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2717 M = header[1]; N = header[2]; nz = header[3]; 2718 2719 if (nz < 0) { 2720 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2721 } 2722 2723 /* read in row lengths */ 2724 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2725 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2726 2727 /* create our matrix */ 2728 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A);CHKERRQ(ierr); 2729 B = *A; 2730 a = (Mat_SeqAIJ*)B->data; 2731 shift = a->indexshift; 2732 2733 /* read in column indices and adjust for Fortran indexing*/ 2734 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2735 if (shift) { 2736 for (i=0; i<nz; i++) { 2737 a->j[i] += 1; 2738 } 2739 } 2740 2741 /* read in nonzero values */ 2742 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2743 2744 /* set matrix "i" values */ 2745 a->i[0] = -shift; 2746 for (i=1; i<= M; i++) { 2747 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2748 a->ilen[i-1] = rowlengths[i-1]; 2749 } 2750 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2751 2752 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2753 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2754 PetscFunctionReturn(0); 2755 } 2756 EXTERN_C_END 2757 2758 #undef __FUNCT__ 2759 #define __FUNCT__ "MatEqual_SeqAIJ" 2760 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2761 { 2762 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2763 int ierr; 2764 PetscTruth flag; 2765 2766 PetscFunctionBegin; 2767 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flag);CHKERRQ(ierr); 2768 if (!flag) SETERRQ(PETSC_ERR_ARG_INCOMP,"Matrices must be same type"); 2769 2770 /* If the matrix dimensions are not equal,or no of nonzeros or shift */ 2771 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)|| (a->indexshift != b->indexshift)) { 2772 *flg = PETSC_FALSE; 2773 PetscFunctionReturn(0); 2774 } 2775 2776 /* if the a->i are the same */ 2777 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2778 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2779 2780 /* if a->j are the same */ 2781 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2782 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2783 2784 /* if a->a are the same */ 2785 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(Scalar),flg);CHKERRQ(ierr); 2786 2787 PetscFunctionReturn(0); 2788 2789 } 2790 2791 #undef __FUNCT__ 2792 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2793 /*@C 2794 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2795 provided by the user. 2796 2797 Coolective on MPI_Comm 2798 2799 Input Parameters: 2800 + comm - must be an MPI communicator of size 1 2801 . m - number of rows 2802 . n - number of columns 2803 . i - row indices 2804 . j - column indices 2805 - a - matrix values 2806 2807 Output Parameter: 2808 . mat - the matrix 2809 2810 Level: intermediate 2811 2812 Notes: 2813 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2814 once the matrix is destroyed 2815 2816 You cannot set new nonzero locations into this matrix, that will generate an error. 2817 2818 The i and j indices can be either 0- or 1 based 2819 2820 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2821 2822 @*/ 2823 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,Scalar *a,Mat *mat) 2824 { 2825 int ierr,ii; 2826 Mat_SeqAIJ *aij; 2827 2828 PetscFunctionBegin; 2829 ierr = MatCreateSeqAIJ(comm,m,n,0,0,mat);CHKERRQ(ierr); 2830 aij = (Mat_SeqAIJ*)(*mat)->data; 2831 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2832 2833 if (i[0] == 1) { 2834 aij->indexshift = -1; 2835 } else if (i[0]) { 2836 SETERRQ(1,"i (row indices) do not start with 0 or 1"); 2837 } 2838 aij->i = i; 2839 aij->j = j; 2840 aij->a = a; 2841 aij->singlemalloc = PETSC_FALSE; 2842 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2843 aij->freedata = PETSC_FALSE; 2844 2845 for (ii=0; ii<m; ii++) { 2846 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2847 #if defined(PETSC_USE_BOPT_g) 2848 if (i[ii+1] - i[ii] < 0) SETERRQ2(1,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 2849 #endif 2850 } 2851 #if defined(PETSC_USE_BOPT_g) 2852 for (ii=0; ii<aij->i[m]; ii++) { 2853 if (j[ii] < -aij->indexshift) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 2854 if (j[ii] > n - 1 -aij->indexshift) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 2855 } 2856 #endif 2857 2858 /* changes indices to start at 0 */ 2859 if (i[0]) { 2860 aij->indexshift = 0; 2861 for (ii=0; ii<m; ii++) { 2862 i[ii]--; 2863 } 2864 for (ii=0; ii<i[m]; ii++) { 2865 j[ii]--; 2866 } 2867 } 2868 2869 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2870 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2871 PetscFunctionReturn(0); 2872 } 2873 2874 #undef __FUNCT__ 2875 #define __FUNCT__ "MatSetColoring_SeqAIJ" 2876 int MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 2877 { 2878 int ierr; 2879 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2880 2881 PetscFunctionBegin; 2882 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 2883 a->coloring = coloring; 2884 PetscFunctionReturn(0); 2885 } 2886 2887 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) 2888 EXTERN_C_BEGIN 2889 #include "adic_utils.h" 2890 EXTERN_C_END 2891 2892 #undef __FUNCT__ 2893 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 2894 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 2895 { 2896 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2897 int m = A->m,*ii = a->i,*jj = a->j,nz,i,*color,j,nlen; 2898 Scalar *v = a->a,*values; 2899 char *cadvalues = (char *)advalues; 2900 2901 PetscFunctionBegin; 2902 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 2903 nlen = my_AD_GetDerivTypeSize(); 2904 color = a->coloring->colors; 2905 /* loop over rows */ 2906 for (i=0; i<m; i++) { 2907 nz = ii[i+1] - ii[i]; 2908 /* loop over columns putting computed value into matrix */ 2909 values = my_AD_GetGradArray(cadvalues); 2910 for (j=0; j<nz; j++) { 2911 *v++ = values[color[*jj++]]; 2912 } 2913 cadvalues += nlen; /* jump to next row of derivatives */ 2914 } 2915 PetscFunctionReturn(0); 2916 } 2917 2918 #else 2919 2920 #undef __FUNCT__ 2921 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 2922 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 2923 { 2924 PetscFunctionBegin; 2925 SETERRQ(1,"PETSc installed without ADIC"); 2926 } 2927 2928 #endif 2929 2930 #undef __FUNCT__ 2931 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 2932 int MatSetValuesAdifor_SeqAIJ(Mat A,int nl,void *advalues) 2933 { 2934 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2935 int m = A->m,*ii = a->i,*jj = a->j,nz,i,*color,j; 2936 Scalar *v = a->a,*values = (Scalar *)advalues; 2937 2938 PetscFunctionBegin; 2939 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 2940 color = a->coloring->colors; 2941 /* loop over rows */ 2942 for (i=0; i<m; i++) { 2943 nz = ii[i+1] - ii[i]; 2944 /* loop over columns putting computed value into matrix */ 2945 for (j=0; j<nz; j++) { 2946 *v++ = values[color[*jj++]]; 2947 } 2948 values += nl; /* jump to next row of derivatives */ 2949 } 2950 PetscFunctionReturn(0); 2951 } 2952 2953