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