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