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