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