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,isascii,isbinary,isdraw; 562 563 PetscFunctionBegin; 564 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 565 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);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 (isascii) { 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,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 SETERRQ(1,"Incompatible A/B sizes for symmetry test"); 1410 aii = aij->i; bii = bij->i; 1411 adx = aij->j; bdx = bij->j; 1412 va = aij->a; vb = bij->a; 1413 ierr = PetscMalloc(ma*sizeof(int),&aptr);CHKERRQ(ierr); 1414 ierr = PetscMalloc(mb*sizeof(int),&bptr);CHKERRQ(ierr); 1415 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1416 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1417 1418 *f = PETSC_TRUE; 1419 for (i=0; i<ma; i++) { 1420 /*printf("row %d spans %d--%d; we start @ %d\n", 1421 i,idx[ii[i]],idx[ii[i+1]-1],idx[aptr[i]]);*/ 1422 while (aptr[i]<aii[i+1]) { 1423 int idc,idr; PetscScalar vc,vr; 1424 /* column/row index/value */ 1425 idc = adx[aptr[i]]; idr = bdx[bptr[idc]]; 1426 vc = va[aptr[i]]; vr = vb[bptr[idc]]; 1427 /*printf("comparing %d: (%d,%d)=%e to (%d,%d)=%e\n", 1428 aptr[i],i,idc,vc,idc,idr,vr);*/ 1429 if (i!=idr || vc!=vr) { 1430 *f = PETSC_FALSE; goto done; 1431 } else { 1432 aptr[i]++; if (B || i!=idc) bptr[idc]++; 1433 } 1434 } 1435 } 1436 done: 1437 ierr = PetscFree(aptr);CHKERRQ(ierr); 1438 if (B) { 1439 ierr = PetscFree(bptr);CHKERRQ(ierr); 1440 } 1441 1442 PetscFunctionReturn(0); 1443 } 1444 EXTERN_C_END 1445 1446 #undef __FUNCT__ 1447 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1448 int MatIsSymmetric_SeqAIJ(Mat A,PetscTruth *f) 1449 { 1450 int ierr; 1451 PetscFunctionBegin; 1452 ierr = MatIsTranspose_SeqAIJ(A,A,f);CHKERRQ(ierr); 1453 PetscFunctionReturn(0); 1454 } 1455 1456 #undef __FUNCT__ 1457 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1458 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1459 { 1460 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1461 PetscScalar *l,*r,x,*v; 1462 int ierr,i,j,m = A->m,n = A->n,M,nz = a->nz,*jj; 1463 1464 PetscFunctionBegin; 1465 if (ll) { 1466 /* The local size is used so that VecMPI can be passed to this routine 1467 by MatDiagonalScale_MPIAIJ */ 1468 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1469 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1470 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1471 v = a->a; 1472 for (i=0; i<m; i++) { 1473 x = l[i]; 1474 M = a->i[i+1] - a->i[i]; 1475 for (j=0; j<M; j++) { (*v++) *= x;} 1476 } 1477 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1478 PetscLogFlops(nz); 1479 } 1480 if (rr) { 1481 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1482 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1483 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1484 v = a->a; jj = a->j; 1485 for (i=0; i<nz; i++) { 1486 (*v++) *= r[*jj++]; 1487 } 1488 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1489 PetscLogFlops(nz); 1490 } 1491 PetscFunctionReturn(0); 1492 } 1493 1494 #undef __FUNCT__ 1495 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1496 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatReuse scall,Mat *B) 1497 { 1498 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1499 int *smap,i,k,kstart,kend,ierr,oldcols = A->n,*lens; 1500 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1501 int *irow,*icol,nrows,ncols; 1502 int *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1503 PetscScalar *a_new,*mat_a; 1504 Mat C; 1505 PetscTruth stride; 1506 1507 PetscFunctionBegin; 1508 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1509 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1510 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1511 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1512 1513 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1514 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1515 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1516 1517 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1518 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1519 if (stride && step == 1) { 1520 /* special case of contiguous rows */ 1521 ierr = PetscMalloc((2*nrows+1)*sizeof(int),&lens);CHKERRQ(ierr); 1522 starts = lens + nrows; 1523 /* loop over new rows determining lens and starting points */ 1524 for (i=0; i<nrows; i++) { 1525 kstart = ai[irow[i]]; 1526 kend = kstart + ailen[irow[i]]; 1527 for (k=kstart; k<kend; k++) { 1528 if (aj[k] >= first) { 1529 starts[i] = k; 1530 break; 1531 } 1532 } 1533 sum = 0; 1534 while (k < kend) { 1535 if (aj[k++] >= first+ncols) break; 1536 sum++; 1537 } 1538 lens[i] = sum; 1539 } 1540 /* create submatrix */ 1541 if (scall == MAT_REUSE_MATRIX) { 1542 int n_cols,n_rows; 1543 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1544 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1545 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1546 C = *B; 1547 } else { 1548 ierr = MatCreate(A->comm,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE,&C);CHKERRQ(ierr); 1549 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1550 ierr = MatSeqAIJSetPreallocation(C,0,lens);CHKERRQ(ierr); 1551 } 1552 c = (Mat_SeqAIJ*)C->data; 1553 1554 /* loop over rows inserting into submatrix */ 1555 a_new = c->a; 1556 j_new = c->j; 1557 i_new = c->i; 1558 1559 for (i=0; i<nrows; i++) { 1560 ii = starts[i]; 1561 lensi = lens[i]; 1562 for (k=0; k<lensi; k++) { 1563 *j_new++ = aj[ii+k] - first; 1564 } 1565 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1566 a_new += lensi; 1567 i_new[i+1] = i_new[i] + lensi; 1568 c->ilen[i] = lensi; 1569 } 1570 ierr = PetscFree(lens);CHKERRQ(ierr); 1571 } else { 1572 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1573 ierr = PetscMalloc((1+oldcols)*sizeof(int),&smap);CHKERRQ(ierr); 1574 1575 ierr = PetscMalloc((1+nrows)*sizeof(int),&lens);CHKERRQ(ierr); 1576 ierr = PetscMemzero(smap,oldcols*sizeof(int));CHKERRQ(ierr); 1577 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1578 /* determine lens of each row */ 1579 for (i=0; i<nrows; i++) { 1580 kstart = ai[irow[i]]; 1581 kend = kstart + a->ilen[irow[i]]; 1582 lens[i] = 0; 1583 for (k=kstart; k<kend; k++) { 1584 if (smap[aj[k]]) { 1585 lens[i]++; 1586 } 1587 } 1588 } 1589 /* Create and fill new matrix */ 1590 if (scall == MAT_REUSE_MATRIX) { 1591 PetscTruth equal; 1592 1593 c = (Mat_SeqAIJ *)((*B)->data); 1594 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1595 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(int),&equal);CHKERRQ(ierr); 1596 if (!equal) { 1597 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1598 } 1599 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(int));CHKERRQ(ierr); 1600 C = *B; 1601 } else { 1602 ierr = MatCreate(A->comm,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE,&C);CHKERRQ(ierr); 1603 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1604 ierr = MatSeqAIJSetPreallocation(C,0,lens);CHKERRQ(ierr); 1605 } 1606 c = (Mat_SeqAIJ *)(C->data); 1607 for (i=0; i<nrows; i++) { 1608 row = irow[i]; 1609 kstart = ai[row]; 1610 kend = kstart + a->ilen[row]; 1611 mat_i = c->i[i]; 1612 mat_j = c->j + mat_i; 1613 mat_a = c->a + mat_i; 1614 mat_ilen = c->ilen + i; 1615 for (k=kstart; k<kend; k++) { 1616 if ((tcol=smap[a->j[k]])) { 1617 *mat_j++ = tcol - 1; 1618 *mat_a++ = a->a[k]; 1619 (*mat_ilen)++; 1620 1621 } 1622 } 1623 } 1624 /* Free work space */ 1625 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1626 ierr = PetscFree(smap);CHKERRQ(ierr); 1627 ierr = PetscFree(lens);CHKERRQ(ierr); 1628 } 1629 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1630 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1631 1632 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1633 *B = C; 1634 PetscFunctionReturn(0); 1635 } 1636 1637 /* 1638 */ 1639 #undef __FUNCT__ 1640 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1641 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1642 { 1643 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1644 int ierr; 1645 Mat outA; 1646 PetscTruth row_identity,col_identity; 1647 1648 PetscFunctionBegin; 1649 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1650 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1651 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1652 if (!row_identity || !col_identity) { 1653 SETERRQ(1,"Row and column permutations must be identity for in-place ILU"); 1654 } 1655 1656 outA = inA; 1657 inA->factor = FACTOR_LU; 1658 a->row = row; 1659 a->col = col; 1660 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1661 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1662 1663 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1664 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1665 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1666 PetscLogObjectParent(inA,a->icol); 1667 1668 if (!a->solve_work) { /* this matrix may have been factored before */ 1669 ierr = PetscMalloc((inA->m+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1670 } 1671 1672 if (!a->diag) { 1673 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1674 } 1675 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1676 PetscFunctionReturn(0); 1677 } 1678 1679 #include "petscblaslapack.h" 1680 #undef __FUNCT__ 1681 #define __FUNCT__ "MatScale_SeqAIJ" 1682 int MatScale_SeqAIJ(const PetscScalar *alpha,Mat inA) 1683 { 1684 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1685 int one = 1; 1686 1687 PetscFunctionBegin; 1688 BLscal_(&a->nz,(PetscScalar*)alpha,a->a,&one); 1689 PetscLogFlops(a->nz); 1690 PetscFunctionReturn(0); 1691 } 1692 1693 #undef __FUNCT__ 1694 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1695 int MatGetSubMatrices_SeqAIJ(Mat A,int n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1696 { 1697 int ierr,i; 1698 1699 PetscFunctionBegin; 1700 if (scall == MAT_INITIAL_MATRIX) { 1701 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1702 } 1703 1704 for (i=0; i<n; i++) { 1705 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1706 } 1707 PetscFunctionReturn(0); 1708 } 1709 1710 #undef __FUNCT__ 1711 #define __FUNCT__ "MatGetBlockSize_SeqAIJ" 1712 int MatGetBlockSize_SeqAIJ(Mat A,int *bs) 1713 { 1714 PetscFunctionBegin; 1715 *bs = 1; 1716 PetscFunctionReturn(0); 1717 } 1718 1719 #undef __FUNCT__ 1720 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1721 int MatIncreaseOverlap_SeqAIJ(Mat A,int is_max,IS is[],int ov) 1722 { 1723 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1724 int row,i,j,k,l,m,n,*idx,ierr,*nidx,isz,val; 1725 int start,end,*ai,*aj; 1726 PetscBT table; 1727 1728 PetscFunctionBegin; 1729 m = A->m; 1730 ai = a->i; 1731 aj = a->j; 1732 1733 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1734 1735 ierr = PetscMalloc((m+1)*sizeof(int),&nidx);CHKERRQ(ierr); 1736 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1737 1738 for (i=0; i<is_max; i++) { 1739 /* Initialize the two local arrays */ 1740 isz = 0; 1741 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1742 1743 /* Extract the indices, assume there can be duplicate entries */ 1744 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1745 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1746 1747 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1748 for (j=0; j<n ; ++j){ 1749 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1750 } 1751 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1752 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1753 1754 k = 0; 1755 for (j=0; j<ov; j++){ /* for each overlap */ 1756 n = isz; 1757 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1758 row = nidx[k]; 1759 start = ai[row]; 1760 end = ai[row+1]; 1761 for (l = start; l<end ; l++){ 1762 val = aj[l] ; 1763 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1764 } 1765 } 1766 } 1767 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1768 } 1769 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1770 ierr = PetscFree(nidx);CHKERRQ(ierr); 1771 PetscFunctionReturn(0); 1772 } 1773 1774 /* -------------------------------------------------------------- */ 1775 #undef __FUNCT__ 1776 #define __FUNCT__ "MatPermute_SeqAIJ" 1777 int MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1778 { 1779 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1780 PetscScalar *vwork; 1781 int i,ierr,nz,m = A->m,n = A->n,*cwork; 1782 int *row,*col,*cnew,j,*lens; 1783 IS icolp,irowp; 1784 1785 PetscFunctionBegin; 1786 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1787 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1788 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1789 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1790 1791 /* determine lengths of permuted rows */ 1792 ierr = PetscMalloc((m+1)*sizeof(int),&lens);CHKERRQ(ierr); 1793 for (i=0; i<m; i++) { 1794 lens[row[i]] = a->i[i+1] - a->i[i]; 1795 } 1796 ierr = MatCreate(A->comm,m,n,m,n,B);CHKERRQ(ierr); 1797 ierr = MatSetType(*B,A->type_name);CHKERRQ(ierr); 1798 ierr = MatSeqAIJSetPreallocation(*B,0,lens);CHKERRQ(ierr); 1799 ierr = PetscFree(lens);CHKERRQ(ierr); 1800 1801 ierr = PetscMalloc(n*sizeof(int),&cnew);CHKERRQ(ierr); 1802 for (i=0; i<m; i++) { 1803 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1804 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1805 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1806 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1807 } 1808 ierr = PetscFree(cnew);CHKERRQ(ierr); 1809 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1810 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1811 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1812 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1813 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1814 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1815 PetscFunctionReturn(0); 1816 } 1817 1818 #undef __FUNCT__ 1819 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1820 int MatPrintHelp_SeqAIJ(Mat A) 1821 { 1822 static PetscTruth called = PETSC_FALSE; 1823 MPI_Comm comm = A->comm; 1824 int ierr; 1825 1826 PetscFunctionBegin; 1827 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1828 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1829 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1830 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1831 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1832 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1833 PetscFunctionReturn(0); 1834 } 1835 1836 #undef __FUNCT__ 1837 #define __FUNCT__ "MatCopy_SeqAIJ" 1838 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1839 { 1840 int ierr; 1841 1842 PetscFunctionBegin; 1843 /* If the two matrices have the same copy implementation, use fast copy. */ 1844 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1845 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1846 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1847 1848 if (a->i[A->m] != b->i[B->m]) { 1849 SETERRQ(1,"Number of nonzeros in two matrices are different"); 1850 } 1851 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 1852 } else { 1853 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1854 } 1855 PetscFunctionReturn(0); 1856 } 1857 1858 #undef __FUNCT__ 1859 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1860 int MatSetUpPreallocation_SeqAIJ(Mat A) 1861 { 1862 int ierr; 1863 1864 PetscFunctionBegin; 1865 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1866 PetscFunctionReturn(0); 1867 } 1868 1869 #undef __FUNCT__ 1870 #define __FUNCT__ "MatGetArray_SeqAIJ" 1871 int MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 1872 { 1873 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1874 PetscFunctionBegin; 1875 *array = a->a; 1876 PetscFunctionReturn(0); 1877 } 1878 1879 #undef __FUNCT__ 1880 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1881 int MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 1882 { 1883 PetscFunctionBegin; 1884 PetscFunctionReturn(0); 1885 } 1886 1887 #undef __FUNCT__ 1888 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1889 int MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1890 { 1891 int (*f)(void *,Vec,Vec,void*) = (int (*)(void *,Vec,Vec,void *))coloring->f; 1892 int k,ierr,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1893 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1894 PetscScalar *vscale_array; 1895 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1896 Vec w1,w2,w3; 1897 void *fctx = coloring->fctx; 1898 PetscTruth flg; 1899 1900 PetscFunctionBegin; 1901 if (!coloring->w1) { 1902 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1903 PetscLogObjectParent(coloring,coloring->w1); 1904 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1905 PetscLogObjectParent(coloring,coloring->w2); 1906 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1907 PetscLogObjectParent(coloring,coloring->w3); 1908 } 1909 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1910 1911 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1912 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1913 if (flg) { 1914 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 1915 } else { 1916 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1917 } 1918 1919 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1920 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1921 1922 /* 1923 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1924 coloring->F for the coarser grids from the finest 1925 */ 1926 if (coloring->F) { 1927 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1928 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1929 if (m1 != m2) { 1930 coloring->F = 0; 1931 } 1932 } 1933 1934 if (coloring->F) { 1935 w1 = coloring->F; 1936 coloring->F = 0; 1937 } else { 1938 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1939 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1940 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1941 } 1942 1943 /* 1944 Compute all the scale factors and share with other processors 1945 */ 1946 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 1947 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 1948 for (k=0; k<coloring->ncolors; k++) { 1949 /* 1950 Loop over each column associated with color adding the 1951 perturbation to the vector w3. 1952 */ 1953 for (l=0; l<coloring->ncolumns[k]; l++) { 1954 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1955 dx = xx[col]; 1956 if (dx == 0.0) dx = 1.0; 1957 #if !defined(PETSC_USE_COMPLEX) 1958 if (dx < umin && dx >= 0.0) dx = umin; 1959 else if (dx < 0.0 && dx > -umin) dx = -umin; 1960 #else 1961 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1962 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1963 #endif 1964 dx *= epsilon; 1965 vscale_array[col] = 1.0/dx; 1966 } 1967 } 1968 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1969 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1970 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1971 1972 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 1973 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 1974 1975 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 1976 else vscaleforrow = coloring->columnsforrow; 1977 1978 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1979 /* 1980 Loop over each color 1981 */ 1982 for (k=0; k<coloring->ncolors; k++) { 1983 coloring->currentcolor = k; 1984 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 1985 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 1986 /* 1987 Loop over each column associated with color adding the 1988 perturbation to the vector w3. 1989 */ 1990 for (l=0; l<coloring->ncolumns[k]; l++) { 1991 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1992 dx = xx[col]; 1993 if (dx == 0.0) dx = 1.0; 1994 #if !defined(PETSC_USE_COMPLEX) 1995 if (dx < umin && dx >= 0.0) dx = umin; 1996 else if (dx < 0.0 && dx > -umin) dx = -umin; 1997 #else 1998 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1999 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2000 #endif 2001 dx *= epsilon; 2002 if (!PetscAbsScalar(dx)) SETERRQ(1,"Computed 0 differencing parameter"); 2003 w3_array[col] += dx; 2004 } 2005 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2006 2007 /* 2008 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2009 */ 2010 2011 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2012 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2013 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2014 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 2015 2016 /* 2017 Loop over rows of vector, putting results into Jacobian matrix 2018 */ 2019 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2020 for (l=0; l<coloring->nrows[k]; l++) { 2021 row = coloring->rows[k][l]; 2022 col = coloring->columnsforrow[k][l]; 2023 y[row] *= vscale_array[vscaleforrow[k][l]]; 2024 srow = row + start; 2025 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2026 } 2027 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2028 } 2029 coloring->currentcolor = k; 2030 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2031 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2032 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2033 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2034 PetscFunctionReturn(0); 2035 } 2036 2037 #include "petscblaslapack.h" 2038 #undef __FUNCT__ 2039 #define __FUNCT__ "MatAXPY_SeqAIJ" 2040 int MatAXPY_SeqAIJ(const PetscScalar a[],Mat X,Mat Y,MatStructure str) 2041 { 2042 int ierr,one=1,i; 2043 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2044 2045 PetscFunctionBegin; 2046 if (str == SAME_NONZERO_PATTERN) { 2047 BLaxpy_(&x->nz,(PetscScalar*)a,x->a,&one,y->a,&one); 2048 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2049 if (y->xtoy && y->XtoY != X) { 2050 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2051 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2052 } 2053 if (!y->xtoy) { /* get xtoy */ 2054 ierr = MatAXPYGetxtoy_Private(X->m,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2055 y->XtoY = X; 2056 } 2057 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += (*a)*(x->a[i]); 2058 PetscLogInfo(0,"MatAXPY_SeqAIJ: ratio of nnz(X)/nnz(Y): %d/%d = %g\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz); 2059 } else { 2060 ierr = MatAXPY_Basic(a,X,Y,str);CHKERRQ(ierr); 2061 } 2062 PetscFunctionReturn(0); 2063 } 2064 2065 /* -------------------------------------------------------------------*/ 2066 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2067 MatGetRow_SeqAIJ, 2068 MatRestoreRow_SeqAIJ, 2069 MatMult_SeqAIJ, 2070 /* 4*/ MatMultAdd_SeqAIJ, 2071 MatMultTranspose_SeqAIJ, 2072 MatMultTransposeAdd_SeqAIJ, 2073 MatSolve_SeqAIJ, 2074 MatSolveAdd_SeqAIJ, 2075 MatSolveTranspose_SeqAIJ, 2076 /*10*/ MatSolveTransposeAdd_SeqAIJ, 2077 MatLUFactor_SeqAIJ, 2078 0, 2079 MatRelax_SeqAIJ, 2080 MatTranspose_SeqAIJ, 2081 /*15*/ MatGetInfo_SeqAIJ, 2082 MatEqual_SeqAIJ, 2083 MatGetDiagonal_SeqAIJ, 2084 MatDiagonalScale_SeqAIJ, 2085 MatNorm_SeqAIJ, 2086 /*20*/ 0, 2087 MatAssemblyEnd_SeqAIJ, 2088 MatCompress_SeqAIJ, 2089 MatSetOption_SeqAIJ, 2090 MatZeroEntries_SeqAIJ, 2091 /*25*/ MatZeroRows_SeqAIJ, 2092 MatLUFactorSymbolic_SeqAIJ, 2093 MatLUFactorNumeric_SeqAIJ, 2094 MatCholeskyFactorSymbolic_SeqAIJ, 2095 MatCholeskyFactorNumeric_SeqAIJ, 2096 /*30*/ MatSetUpPreallocation_SeqAIJ, 2097 MatILUFactorSymbolic_SeqAIJ, 2098 MatICCFactorSymbolic_SeqAIJ, 2099 MatGetArray_SeqAIJ, 2100 MatRestoreArray_SeqAIJ, 2101 /*35*/ MatDuplicate_SeqAIJ, 2102 0, 2103 0, 2104 MatILUFactor_SeqAIJ, 2105 0, 2106 /*40*/ MatAXPY_SeqAIJ, 2107 MatGetSubMatrices_SeqAIJ, 2108 MatIncreaseOverlap_SeqAIJ, 2109 MatGetValues_SeqAIJ, 2110 MatCopy_SeqAIJ, 2111 /*45*/ MatPrintHelp_SeqAIJ, 2112 MatScale_SeqAIJ, 2113 0, 2114 0, 2115 MatILUDTFactor_SeqAIJ, 2116 /*50*/ MatGetBlockSize_SeqAIJ, 2117 MatGetRowIJ_SeqAIJ, 2118 MatRestoreRowIJ_SeqAIJ, 2119 MatGetColumnIJ_SeqAIJ, 2120 MatRestoreColumnIJ_SeqAIJ, 2121 /*55*/ MatFDColoringCreate_SeqAIJ, 2122 0, 2123 0, 2124 MatPermute_SeqAIJ, 2125 0, 2126 /*60*/ 0, 2127 MatDestroy_SeqAIJ, 2128 MatView_SeqAIJ, 2129 MatGetPetscMaps_Petsc, 2130 0, 2131 /*65*/ 0, 2132 0, 2133 0, 2134 0, 2135 0, 2136 /*70*/ 0, 2137 0, 2138 MatSetColoring_SeqAIJ, 2139 MatSetValuesAdic_SeqAIJ, 2140 MatSetValuesAdifor_SeqAIJ, 2141 /*75*/ MatFDColoringApply_SeqAIJ, 2142 0, 2143 0, 2144 0, 2145 0, 2146 /*80*/ 0, 2147 0, 2148 0, 2149 0, 2150 /*85*/ MatLoad_SeqAIJ, 2151 MatIsSymmetric_SeqAIJ, 2152 }; 2153 2154 EXTERN_C_BEGIN 2155 #undef __FUNCT__ 2156 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2157 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 2158 { 2159 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2160 int i,nz,n; 2161 2162 PetscFunctionBegin; 2163 2164 nz = aij->maxnz; 2165 n = mat->n; 2166 for (i=0; i<nz; i++) { 2167 aij->j[i] = indices[i]; 2168 } 2169 aij->nz = nz; 2170 for (i=0; i<n; i++) { 2171 aij->ilen[i] = aij->imax[i]; 2172 } 2173 2174 PetscFunctionReturn(0); 2175 } 2176 EXTERN_C_END 2177 2178 #undef __FUNCT__ 2179 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2180 /*@ 2181 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2182 in the matrix. 2183 2184 Input Parameters: 2185 + mat - the SeqAIJ matrix 2186 - indices - the column indices 2187 2188 Level: advanced 2189 2190 Notes: 2191 This can be called if you have precomputed the nonzero structure of the 2192 matrix and want to provide it to the matrix object to improve the performance 2193 of the MatSetValues() operation. 2194 2195 You MUST have set the correct numbers of nonzeros per row in the call to 2196 MatCreateSeqAIJ(). 2197 2198 MUST be called before any calls to MatSetValues(); 2199 2200 The indices should start with zero, not one. 2201 2202 @*/ 2203 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 2204 { 2205 int ierr,(*f)(Mat,int *); 2206 2207 PetscFunctionBegin; 2208 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2209 PetscValidPointer(indices,2); 2210 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2211 if (f) { 2212 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2213 } else { 2214 SETERRQ(1,"Wrong type of matrix to set column indices"); 2215 } 2216 PetscFunctionReturn(0); 2217 } 2218 2219 /* ----------------------------------------------------------------------------------------*/ 2220 2221 EXTERN_C_BEGIN 2222 #undef __FUNCT__ 2223 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2224 int MatStoreValues_SeqAIJ(Mat mat) 2225 { 2226 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2227 size_t nz = aij->i[mat->m],ierr; 2228 2229 PetscFunctionBegin; 2230 if (aij->nonew != 1) { 2231 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2232 } 2233 2234 /* allocate space for values if not already there */ 2235 if (!aij->saved_values) { 2236 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2237 } 2238 2239 /* copy values over */ 2240 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2241 PetscFunctionReturn(0); 2242 } 2243 EXTERN_C_END 2244 2245 #undef __FUNCT__ 2246 #define __FUNCT__ "MatStoreValues" 2247 /*@ 2248 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2249 example, reuse of the linear part of a Jacobian, while recomputing the 2250 nonlinear portion. 2251 2252 Collect on Mat 2253 2254 Input Parameters: 2255 . mat - the matrix (currently on AIJ matrices support this option) 2256 2257 Level: advanced 2258 2259 Common Usage, with SNESSolve(): 2260 $ Create Jacobian matrix 2261 $ Set linear terms into matrix 2262 $ Apply boundary conditions to matrix, at this time matrix must have 2263 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2264 $ boundary conditions again will not change the nonzero structure 2265 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2266 $ ierr = MatStoreValues(mat); 2267 $ Call SNESSetJacobian() with matrix 2268 $ In your Jacobian routine 2269 $ ierr = MatRetrieveValues(mat); 2270 $ Set nonlinear terms in matrix 2271 2272 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2273 $ // build linear portion of Jacobian 2274 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2275 $ ierr = MatStoreValues(mat); 2276 $ loop over nonlinear iterations 2277 $ ierr = MatRetrieveValues(mat); 2278 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2279 $ // call MatAssemblyBegin/End() on matrix 2280 $ Solve linear system with Jacobian 2281 $ endloop 2282 2283 Notes: 2284 Matrix must already be assemblied before calling this routine 2285 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2286 calling this routine. 2287 2288 .seealso: MatRetrieveValues() 2289 2290 @*/ 2291 int MatStoreValues(Mat mat) 2292 { 2293 int ierr,(*f)(Mat); 2294 2295 PetscFunctionBegin; 2296 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2297 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2298 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2299 2300 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2301 if (f) { 2302 ierr = (*f)(mat);CHKERRQ(ierr); 2303 } else { 2304 SETERRQ(1,"Wrong type of matrix to store values"); 2305 } 2306 PetscFunctionReturn(0); 2307 } 2308 2309 EXTERN_C_BEGIN 2310 #undef __FUNCT__ 2311 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2312 int MatRetrieveValues_SeqAIJ(Mat mat) 2313 { 2314 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2315 int nz = aij->i[mat->m],ierr; 2316 2317 PetscFunctionBegin; 2318 if (aij->nonew != 1) { 2319 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2320 } 2321 if (!aij->saved_values) { 2322 SETERRQ(1,"Must call MatStoreValues(A);first"); 2323 } 2324 2325 /* copy values over */ 2326 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2327 PetscFunctionReturn(0); 2328 } 2329 EXTERN_C_END 2330 2331 #undef __FUNCT__ 2332 #define __FUNCT__ "MatRetrieveValues" 2333 /*@ 2334 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2335 example, reuse of the linear part of a Jacobian, while recomputing the 2336 nonlinear portion. 2337 2338 Collect on Mat 2339 2340 Input Parameters: 2341 . mat - the matrix (currently on AIJ matrices support this option) 2342 2343 Level: advanced 2344 2345 .seealso: MatStoreValues() 2346 2347 @*/ 2348 int MatRetrieveValues(Mat mat) 2349 { 2350 int ierr,(*f)(Mat); 2351 2352 PetscFunctionBegin; 2353 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2354 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2355 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2356 2357 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2358 if (f) { 2359 ierr = (*f)(mat);CHKERRQ(ierr); 2360 } else { 2361 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2362 } 2363 PetscFunctionReturn(0); 2364 } 2365 2366 2367 /* --------------------------------------------------------------------------------*/ 2368 #undef __FUNCT__ 2369 #define __FUNCT__ "MatCreateSeqAIJ" 2370 /*@C 2371 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2372 (the default parallel PETSc format). For good matrix assembly performance 2373 the user should preallocate the matrix storage by setting the parameter nz 2374 (or the array nnz). By setting these parameters accurately, performance 2375 during matrix assembly can be increased by more than a factor of 50. 2376 2377 Collective on MPI_Comm 2378 2379 Input Parameters: 2380 + comm - MPI communicator, set to PETSC_COMM_SELF 2381 . m - number of rows 2382 . n - number of columns 2383 . nz - number of nonzeros per row (same for all rows) 2384 - nnz - array containing the number of nonzeros in the various rows 2385 (possibly different for each row) or PETSC_NULL 2386 2387 Output Parameter: 2388 . A - the matrix 2389 2390 Notes: 2391 The AIJ format (also called the Yale sparse matrix format or 2392 compressed row storage), is fully compatible with standard Fortran 77 2393 storage. That is, the stored row and column indices can begin at 2394 either one (as in Fortran) or zero. See the users' manual for details. 2395 2396 Specify the preallocated storage with either nz or nnz (not both). 2397 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2398 allocation. For large problems you MUST preallocate memory or you 2399 will get TERRIBLE performance, see the users' manual chapter on matrices. 2400 2401 By default, this format uses inodes (identical nodes) when possible, to 2402 improve numerical efficiency of matrix-vector products and solves. We 2403 search for consecutive rows with the same nonzero structure, thereby 2404 reusing matrix information to achieve increased efficiency. 2405 2406 Options Database Keys: 2407 + -mat_aij_no_inode - Do not use inodes 2408 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2409 - -mat_aij_oneindex - Internally use indexing starting at 1 2410 rather than 0. Note that when calling MatSetValues(), 2411 the user still MUST index entries starting at 0! 2412 2413 Level: intermediate 2414 2415 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2416 2417 @*/ 2418 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,const int nnz[],Mat *A) 2419 { 2420 int ierr; 2421 2422 PetscFunctionBegin; 2423 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2424 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2425 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2426 PetscFunctionReturn(0); 2427 } 2428 2429 #define SKIP_ALLOCATION -4 2430 2431 #undef __FUNCT__ 2432 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2433 /*@C 2434 MatSeqAIJSetPreallocation - For good matrix assembly performance 2435 the user should preallocate the matrix storage by setting the parameter nz 2436 (or the array nnz). By setting these parameters accurately, performance 2437 during matrix assembly can be increased by more than a factor of 50. 2438 2439 Collective on MPI_Comm 2440 2441 Input Parameters: 2442 + comm - MPI communicator, set to PETSC_COMM_SELF 2443 . m - number of rows 2444 . n - number of columns 2445 . nz - number of nonzeros per row (same for all rows) 2446 - nnz - array containing the number of nonzeros in the various rows 2447 (possibly different for each row) or PETSC_NULL 2448 2449 Output Parameter: 2450 . A - the matrix 2451 2452 Notes: 2453 The AIJ format (also called the Yale sparse matrix format or 2454 compressed row storage), is fully compatible with standard Fortran 77 2455 storage. That is, the stored row and column indices can begin at 2456 either one (as in Fortran) or zero. See the users' manual for details. 2457 2458 Specify the preallocated storage with either nz or nnz (not both). 2459 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2460 allocation. For large problems you MUST preallocate memory or you 2461 will get TERRIBLE performance, see the users' manual chapter on matrices. 2462 2463 By default, this format uses inodes (identical nodes) when possible, to 2464 improve numerical efficiency of matrix-vector products and solves. We 2465 search for consecutive rows with the same nonzero structure, thereby 2466 reusing matrix information to achieve increased efficiency. 2467 2468 Options Database Keys: 2469 + -mat_aij_no_inode - Do not use inodes 2470 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2471 - -mat_aij_oneindex - Internally use indexing starting at 1 2472 rather than 0. Note that when calling MatSetValues(), 2473 the user still MUST index entries starting at 0! 2474 2475 Level: intermediate 2476 2477 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2478 2479 @*/ 2480 int MatSeqAIJSetPreallocation(Mat B,int nz,const int nnz[]) 2481 { 2482 int ierr,(*f)(Mat,int,const int[]); 2483 2484 PetscFunctionBegin; 2485 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2486 if (f) { 2487 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2488 } 2489 PetscFunctionReturn(0); 2490 } 2491 2492 EXTERN_C_BEGIN 2493 #undef __FUNCT__ 2494 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2495 int MatSeqAIJSetPreallocation_SeqAIJ(Mat B,int nz,int *nnz) 2496 { 2497 Mat_SeqAIJ *b; 2498 size_t len = 0; 2499 PetscTruth skipallocation = PETSC_FALSE; 2500 int i,ierr; 2501 2502 PetscFunctionBegin; 2503 2504 if (nz == SKIP_ALLOCATION) { 2505 skipallocation = PETSC_TRUE; 2506 nz = 0; 2507 } 2508 2509 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2510 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2511 if (nnz) { 2512 for (i=0; i<B->m; i++) { 2513 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2514 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); 2515 } 2516 } 2517 2518 B->preallocated = PETSC_TRUE; 2519 b = (Mat_SeqAIJ*)B->data; 2520 2521 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2522 if (!nnz) { 2523 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2524 else if (nz <= 0) nz = 1; 2525 for (i=0; i<B->m; i++) b->imax[i] = nz; 2526 nz = nz*B->m; 2527 } else { 2528 nz = 0; 2529 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2530 } 2531 2532 if (!skipallocation) { 2533 /* allocate the matrix space */ 2534 len = ((size_t) nz)*(sizeof(int) + sizeof(PetscScalar)) + (B->m+1)*sizeof(int); 2535 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2536 b->j = (int*)(b->a + nz); 2537 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2538 b->i = b->j + nz; 2539 b->i[0] = 0; 2540 for (i=1; i<B->m+1; i++) { 2541 b->i[i] = b->i[i-1] + b->imax[i-1]; 2542 } 2543 b->singlemalloc = PETSC_TRUE; 2544 b->freedata = PETSC_TRUE; 2545 } else { 2546 b->freedata = PETSC_FALSE; 2547 } 2548 2549 /* b->ilen will count nonzeros in each row so far. */ 2550 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2551 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2552 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2553 2554 b->nz = 0; 2555 b->maxnz = nz; 2556 B->info.nz_unneeded = (double)b->maxnz; 2557 PetscFunctionReturn(0); 2558 } 2559 EXTERN_C_END 2560 2561 /*MC 2562 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 2563 based on compressed sparse row format. 2564 2565 Options Database Keys: 2566 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 2567 2568 Level: beginner 2569 2570 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 2571 M*/ 2572 2573 EXTERN_C_BEGIN 2574 #undef __FUNCT__ 2575 #define __FUNCT__ "MatCreate_SeqAIJ" 2576 int MatCreate_SeqAIJ(Mat B) 2577 { 2578 Mat_SeqAIJ *b; 2579 int ierr,size; 2580 2581 PetscFunctionBegin; 2582 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2583 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2584 2585 B->m = B->M = PetscMax(B->m,B->M); 2586 B->n = B->N = PetscMax(B->n,B->N); 2587 2588 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2589 B->data = (void*)b; 2590 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2591 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2592 B->factor = 0; 2593 B->lupivotthreshold = 1.0; 2594 B->mapping = 0; 2595 ierr = PetscOptionsGetReal(B->prefix,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2596 ierr = PetscOptionsHasName(B->prefix,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2597 b->row = 0; 2598 b->col = 0; 2599 b->icol = 0; 2600 b->reallocs = 0; 2601 2602 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2603 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2604 2605 b->sorted = PETSC_FALSE; 2606 b->ignorezeroentries = PETSC_FALSE; 2607 b->roworiented = PETSC_TRUE; 2608 b->nonew = 0; 2609 b->diag = 0; 2610 b->solve_work = 0; 2611 B->spptr = 0; 2612 b->inode.use = PETSC_TRUE; 2613 b->inode.node_count = 0; 2614 b->inode.size = 0; 2615 b->inode.limit = 5; 2616 b->inode.max_limit = 5; 2617 b->saved_values = 0; 2618 b->idiag = 0; 2619 b->ssor = 0; 2620 b->keepzeroedrows = PETSC_FALSE; 2621 b->xtoy = 0; 2622 b->XtoY = 0; 2623 2624 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2625 2626 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2627 "MatSeqAIJSetColumnIndices_SeqAIJ", 2628 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2629 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2630 "MatStoreValues_SeqAIJ", 2631 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2632 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2633 "MatRetrieveValues_SeqAIJ", 2634 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2635 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2636 "MatConvert_SeqAIJ_SeqSBAIJ", 2637 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2638 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2639 "MatConvert_SeqAIJ_SeqBAIJ", 2640 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2641 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 2642 "MatIsTranspose_SeqAIJ", 2643 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 2644 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2645 "MatSeqAIJSetPreallocation_SeqAIJ", 2646 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2647 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 2648 "MatReorderForNonzeroDiagonal_SeqAIJ", 2649 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 2650 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatAdjustForInodes_C", 2651 "MatAdjustForInodes_SeqAIJ", 2652 MatAdjustForInodes_SeqAIJ);CHKERRQ(ierr); 2653 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJGetInodeSizes_C", 2654 "MatSeqAIJGetInodeSizes_SeqAIJ", 2655 MatSeqAIJGetInodeSizes_SeqAIJ);CHKERRQ(ierr); 2656 ierr = RegisterApplyPtAPRoutines_Private(B);CHKERRQ(ierr); 2657 PetscFunctionReturn(0); 2658 } 2659 EXTERN_C_END 2660 2661 #undef __FUNCT__ 2662 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2663 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2664 { 2665 Mat C; 2666 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2667 int i,m = A->m,ierr; 2668 size_t len; 2669 2670 PetscFunctionBegin; 2671 *B = 0; 2672 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2673 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 2674 c = (Mat_SeqAIJ*)C->data; 2675 2676 C->factor = A->factor; 2677 c->row = 0; 2678 c->col = 0; 2679 c->icol = 0; 2680 c->keepzeroedrows = a->keepzeroedrows; 2681 C->assembled = PETSC_TRUE; 2682 2683 C->M = A->m; 2684 C->N = A->n; 2685 2686 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2687 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2688 for (i=0; i<m; i++) { 2689 c->imax[i] = a->imax[i]; 2690 c->ilen[i] = a->ilen[i]; 2691 } 2692 2693 /* allocate the matrix space */ 2694 c->singlemalloc = PETSC_TRUE; 2695 len = ((size_t) (m+1))*sizeof(int)+(a->i[m])*(sizeof(PetscScalar)+sizeof(int)); 2696 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2697 c->j = (int*)(c->a + a->i[m] ); 2698 c->i = c->j + a->i[m]; 2699 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2700 if (m > 0) { 2701 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(int));CHKERRQ(ierr); 2702 if (cpvalues == MAT_COPY_VALUES) { 2703 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2704 } else { 2705 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2706 } 2707 } 2708 2709 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2710 c->sorted = a->sorted; 2711 c->roworiented = a->roworiented; 2712 c->nonew = a->nonew; 2713 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2714 c->saved_values = 0; 2715 c->idiag = 0; 2716 c->ssor = 0; 2717 c->ignorezeroentries = a->ignorezeroentries; 2718 c->freedata = PETSC_TRUE; 2719 2720 if (a->diag) { 2721 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2722 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2723 for (i=0; i<m; i++) { 2724 c->diag[i] = a->diag[i]; 2725 } 2726 } else c->diag = 0; 2727 c->inode.use = a->inode.use; 2728 c->inode.limit = a->inode.limit; 2729 c->inode.max_limit = a->inode.max_limit; 2730 if (a->inode.size){ 2731 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2732 c->inode.node_count = a->inode.node_count; 2733 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2734 } else { 2735 c->inode.size = 0; 2736 c->inode.node_count = 0; 2737 } 2738 c->nz = a->nz; 2739 c->maxnz = a->maxnz; 2740 c->solve_work = 0; 2741 C->preallocated = PETSC_TRUE; 2742 2743 *B = C; 2744 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2745 PetscFunctionReturn(0); 2746 } 2747 2748 #undef __FUNCT__ 2749 #define __FUNCT__ "MatLoad_SeqAIJ" 2750 int MatLoad_SeqAIJ(PetscViewer viewer,const MatType type,Mat *A) 2751 { 2752 Mat_SeqAIJ *a; 2753 Mat B; 2754 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N; 2755 MPI_Comm comm; 2756 2757 PetscFunctionBegin; 2758 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2759 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2760 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2761 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2762 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2763 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2764 M = header[1]; N = header[2]; nz = header[3]; 2765 2766 if (nz < 0) { 2767 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2768 } 2769 2770 /* read in row lengths */ 2771 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2772 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2773 2774 /* create our matrix */ 2775 ierr = MatCreate(comm,PETSC_DECIDE,PETSC_DECIDE,M,N,&B);CHKERRQ(ierr); 2776 ierr = MatSetType(B,type);CHKERRQ(ierr); 2777 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 2778 a = (Mat_SeqAIJ*)B->data; 2779 2780 /* read in column indices and adjust for Fortran indexing*/ 2781 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2782 2783 /* read in nonzero values */ 2784 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2785 2786 /* set matrix "i" values */ 2787 a->i[0] = 0; 2788 for (i=1; i<= M; i++) { 2789 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2790 a->ilen[i-1] = rowlengths[i-1]; 2791 } 2792 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2793 2794 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2795 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2796 *A = B; 2797 PetscFunctionReturn(0); 2798 } 2799 2800 #undef __FUNCT__ 2801 #define __FUNCT__ "MatEqual_SeqAIJ" 2802 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2803 { 2804 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2805 int ierr; 2806 2807 PetscFunctionBegin; 2808 /* If the matrix dimensions are not equal,or no of nonzeros */ 2809 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)) { 2810 *flg = PETSC_FALSE; 2811 PetscFunctionReturn(0); 2812 } 2813 2814 /* if the a->i are the same */ 2815 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2816 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2817 2818 /* if a->j are the same */ 2819 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2820 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2821 2822 /* if a->a are the same */ 2823 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2824 2825 PetscFunctionReturn(0); 2826 2827 } 2828 2829 #undef __FUNCT__ 2830 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2831 /*@C 2832 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2833 provided by the user. 2834 2835 Coolective on MPI_Comm 2836 2837 Input Parameters: 2838 + comm - must be an MPI communicator of size 1 2839 . m - number of rows 2840 . n - number of columns 2841 . i - row indices 2842 . j - column indices 2843 - a - matrix values 2844 2845 Output Parameter: 2846 . mat - the matrix 2847 2848 Level: intermediate 2849 2850 Notes: 2851 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2852 once the matrix is destroyed 2853 2854 You cannot set new nonzero locations into this matrix, that will generate an error. 2855 2856 The i and j indices are 0 based 2857 2858 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2859 2860 @*/ 2861 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,PetscScalar *a,Mat *mat) 2862 { 2863 int ierr,ii; 2864 Mat_SeqAIJ *aij; 2865 2866 PetscFunctionBegin; 2867 ierr = MatCreate(comm,m,n,m,n,mat);CHKERRQ(ierr); 2868 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 2869 ierr = MatSeqAIJSetPreallocation(*mat,SKIP_ALLOCATION,0);CHKERRQ(ierr); 2870 aij = (Mat_SeqAIJ*)(*mat)->data; 2871 2872 if (i[0] != 0) { 2873 SETERRQ(1,"i (row indices) must start with 0"); 2874 } 2875 aij->i = i; 2876 aij->j = j; 2877 aij->a = a; 2878 aij->singlemalloc = PETSC_FALSE; 2879 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2880 aij->freedata = PETSC_FALSE; 2881 2882 for (ii=0; ii<m; ii++) { 2883 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2884 #if defined(PETSC_USE_BOPT_g) 2885 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]); 2886 #endif 2887 } 2888 #if defined(PETSC_USE_BOPT_g) 2889 for (ii=0; ii<aij->i[m]; ii++) { 2890 if (j[ii] < 0) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 2891 if (j[ii] > n - 1) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 2892 } 2893 #endif 2894 2895 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2896 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2897 PetscFunctionReturn(0); 2898 } 2899 2900 #undef __FUNCT__ 2901 #define __FUNCT__ "MatSetColoring_SeqAIJ" 2902 int MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 2903 { 2904 int ierr; 2905 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2906 2907 PetscFunctionBegin; 2908 if (coloring->ctype == IS_COLORING_LOCAL) { 2909 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 2910 a->coloring = coloring; 2911 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 2912 int i,*larray; 2913 ISColoring ocoloring; 2914 ISColoringValue *colors; 2915 2916 /* set coloring for diagonal portion */ 2917 ierr = PetscMalloc((A->n+1)*sizeof(int),&larray);CHKERRQ(ierr); 2918 for (i=0; i<A->n; i++) { 2919 larray[i] = i; 2920 } 2921 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 2922 ierr = PetscMalloc((A->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 2923 for (i=0; i<A->n; i++) { 2924 colors[i] = coloring->colors[larray[i]]; 2925 } 2926 ierr = PetscFree(larray);CHKERRQ(ierr); 2927 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 2928 a->coloring = ocoloring; 2929 } 2930 PetscFunctionReturn(0); 2931 } 2932 2933 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2934 EXTERN_C_BEGIN 2935 #include "adic/ad_utils.h" 2936 EXTERN_C_END 2937 2938 #undef __FUNCT__ 2939 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 2940 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 2941 { 2942 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2943 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j,nlen; 2944 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 2945 ISColoringValue *color; 2946 2947 PetscFunctionBegin; 2948 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 2949 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 2950 color = a->coloring->colors; 2951 /* loop over rows */ 2952 for (i=0; i<m; i++) { 2953 nz = ii[i+1] - ii[i]; 2954 /* loop over columns putting computed value into matrix */ 2955 for (j=0; j<nz; j++) { 2956 *v++ = values[color[*jj++]]; 2957 } 2958 values += nlen; /* jump to next row of derivatives */ 2959 } 2960 PetscFunctionReturn(0); 2961 } 2962 2963 #else 2964 2965 #undef __FUNCT__ 2966 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 2967 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 2968 { 2969 PetscFunctionBegin; 2970 SETERRQ(1,"PETSc installed without ADIC"); 2971 } 2972 2973 #endif 2974 2975 #undef __FUNCT__ 2976 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 2977 int MatSetValuesAdifor_SeqAIJ(Mat A,int nl,void *advalues) 2978 { 2979 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2980 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j; 2981 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 2982 ISColoringValue *color; 2983 2984 PetscFunctionBegin; 2985 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 2986 color = a->coloring->colors; 2987 /* loop over rows */ 2988 for (i=0; i<m; i++) { 2989 nz = ii[i+1] - ii[i]; 2990 /* loop over columns putting computed value into matrix */ 2991 for (j=0; j<nz; j++) { 2992 *v++ = values[color[*jj++]]; 2993 } 2994 values += nl; /* jump to next row of derivatives */ 2995 } 2996 PetscFunctionReturn(0); 2997 } 2998 2999