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