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 || format == PETSC_VIEWER_ASCII_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 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 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 ierr = PetscLogFlops(2*a->nz - m);CHKERRQ(ierr); 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 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 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 ierr = PetscLogFlops(m);CHKERRQ(ierr); 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 ierr = PetscLogFlops(2*m);CHKERRQ(ierr); 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 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 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 ierr = PetscLogFlops(6*m-1 + 2*a->nz);CHKERRQ(ierr); 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 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 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 ierr = PetscLogFlops(m);CHKERRQ(ierr); 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 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 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 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 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 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 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 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 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 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 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 PetscErrorCode ierr; 1745 1746 1747 PetscFunctionBegin; 1748 BLASscal_(&bnz,(PetscScalar*)alpha,a->a,&one); 1749 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1750 PetscFunctionReturn(0); 1751 } 1752 1753 #undef __FUNCT__ 1754 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1755 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1756 { 1757 PetscErrorCode ierr; 1758 PetscInt i; 1759 1760 PetscFunctionBegin; 1761 if (scall == MAT_INITIAL_MATRIX) { 1762 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1763 } 1764 1765 for (i=0; i<n; i++) { 1766 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1767 } 1768 PetscFunctionReturn(0); 1769 } 1770 1771 #undef __FUNCT__ 1772 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1773 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1774 { 1775 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1776 PetscErrorCode ierr; 1777 PetscInt row,i,j,k,l,m,n,*idx,*nidx,isz,val; 1778 PetscInt start,end,*ai,*aj; 1779 PetscBT table; 1780 1781 PetscFunctionBegin; 1782 m = A->m; 1783 ai = a->i; 1784 aj = a->j; 1785 1786 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1787 1788 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1789 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1790 1791 for (i=0; i<is_max; i++) { 1792 /* Initialize the two local arrays */ 1793 isz = 0; 1794 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1795 1796 /* Extract the indices, assume there can be duplicate entries */ 1797 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1798 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1799 1800 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1801 for (j=0; j<n ; ++j){ 1802 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1803 } 1804 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1805 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1806 1807 k = 0; 1808 for (j=0; j<ov; j++){ /* for each overlap */ 1809 n = isz; 1810 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1811 row = nidx[k]; 1812 start = ai[row]; 1813 end = ai[row+1]; 1814 for (l = start; l<end ; l++){ 1815 val = aj[l] ; 1816 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1817 } 1818 } 1819 } 1820 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1821 } 1822 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1823 ierr = PetscFree(nidx);CHKERRQ(ierr); 1824 PetscFunctionReturn(0); 1825 } 1826 1827 /* -------------------------------------------------------------- */ 1828 #undef __FUNCT__ 1829 #define __FUNCT__ "MatPermute_SeqAIJ" 1830 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1831 { 1832 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1833 PetscErrorCode ierr; 1834 PetscInt i,nz,m = A->m,n = A->n,*col; 1835 PetscInt *row,*cnew,j,*lens; 1836 IS icolp,irowp; 1837 PetscInt *cwork; 1838 PetscScalar *vwork; 1839 1840 PetscFunctionBegin; 1841 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1842 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1843 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1844 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1845 1846 /* determine lengths of permuted rows */ 1847 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1848 for (i=0; i<m; i++) { 1849 lens[row[i]] = a->i[i+1] - a->i[i]; 1850 } 1851 ierr = MatCreate(A->comm,m,n,m,n,B);CHKERRQ(ierr); 1852 ierr = MatSetType(*B,A->type_name);CHKERRQ(ierr); 1853 ierr = MatSeqAIJSetPreallocation(*B,0,lens);CHKERRQ(ierr); 1854 ierr = PetscFree(lens);CHKERRQ(ierr); 1855 1856 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 1857 for (i=0; i<m; i++) { 1858 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1859 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1860 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1861 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1862 } 1863 ierr = PetscFree(cnew);CHKERRQ(ierr); 1864 (*B)->assembled = PETSC_FALSE; 1865 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1866 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1867 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1868 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1869 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1870 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1871 PetscFunctionReturn(0); 1872 } 1873 1874 #undef __FUNCT__ 1875 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1876 PetscErrorCode MatPrintHelp_SeqAIJ(Mat A) 1877 { 1878 static PetscTruth called = PETSC_FALSE; 1879 MPI_Comm comm = A->comm; 1880 PetscErrorCode ierr; 1881 1882 PetscFunctionBegin; 1883 ierr = MatPrintHelp_Inode(A);CHKERRQ(ierr); 1884 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1885 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1886 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1887 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1888 ierr = (*PetscHelpPrintf)(comm," -mat_no_compressedrow: Do not use compressedrow\n");CHKERRQ(ierr); 1889 PetscFunctionReturn(0); 1890 } 1891 1892 #undef __FUNCT__ 1893 #define __FUNCT__ "MatCopy_SeqAIJ" 1894 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1895 { 1896 PetscErrorCode ierr; 1897 1898 PetscFunctionBegin; 1899 /* If the two matrices have the same copy implementation, use fast copy. */ 1900 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1901 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1902 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1903 1904 if (a->i[A->m] != b->i[B->m]) { 1905 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1906 } 1907 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 1908 } else { 1909 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1910 } 1911 PetscFunctionReturn(0); 1912 } 1913 1914 #undef __FUNCT__ 1915 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1916 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 1917 { 1918 PetscErrorCode ierr; 1919 1920 PetscFunctionBegin; 1921 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1922 PetscFunctionReturn(0); 1923 } 1924 1925 #undef __FUNCT__ 1926 #define __FUNCT__ "MatGetArray_SeqAIJ" 1927 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 1928 { 1929 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1930 PetscFunctionBegin; 1931 *array = a->a; 1932 PetscFunctionReturn(0); 1933 } 1934 1935 #undef __FUNCT__ 1936 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1937 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 1938 { 1939 PetscFunctionBegin; 1940 PetscFunctionReturn(0); 1941 } 1942 1943 #undef __FUNCT__ 1944 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1945 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1946 { 1947 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 1948 PetscErrorCode ierr; 1949 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1950 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1951 PetscScalar *vscale_array; 1952 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1953 Vec w1,w2,w3; 1954 void *fctx = coloring->fctx; 1955 PetscTruth flg; 1956 1957 PetscFunctionBegin; 1958 if (!coloring->w1) { 1959 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1960 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 1961 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1962 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 1963 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1964 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 1965 } 1966 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1967 1968 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1969 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1970 if (flg) { 1971 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 1972 } else { 1973 PetscTruth assembled; 1974 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 1975 if (assembled) { 1976 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1977 } 1978 } 1979 1980 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1981 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1982 1983 /* 1984 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1985 coloring->F for the coarser grids from the finest 1986 */ 1987 if (coloring->F) { 1988 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1989 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1990 if (m1 != m2) { 1991 coloring->F = 0; 1992 } 1993 } 1994 1995 if (coloring->F) { 1996 w1 = coloring->F; 1997 coloring->F = 0; 1998 } else { 1999 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2000 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2001 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2002 } 2003 2004 /* 2005 Compute all the scale factors and share with other processors 2006 */ 2007 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2008 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2009 for (k=0; k<coloring->ncolors; k++) { 2010 /* 2011 Loop over each column associated with color adding the 2012 perturbation to the vector w3. 2013 */ 2014 for (l=0; l<coloring->ncolumns[k]; l++) { 2015 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2016 dx = xx[col]; 2017 if (dx == 0.0) dx = 1.0; 2018 #if !defined(PETSC_USE_COMPLEX) 2019 if (dx < umin && dx >= 0.0) dx = umin; 2020 else if (dx < 0.0 && dx > -umin) dx = -umin; 2021 #else 2022 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2023 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2024 #endif 2025 dx *= epsilon; 2026 vscale_array[col] = 1.0/dx; 2027 } 2028 } 2029 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2030 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2031 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2032 2033 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2034 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2035 2036 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2037 else vscaleforrow = coloring->columnsforrow; 2038 2039 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2040 /* 2041 Loop over each color 2042 */ 2043 for (k=0; k<coloring->ncolors; k++) { 2044 coloring->currentcolor = k; 2045 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2046 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2047 /* 2048 Loop over each column associated with color adding the 2049 perturbation to the vector w3. 2050 */ 2051 for (l=0; l<coloring->ncolumns[k]; l++) { 2052 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2053 dx = xx[col]; 2054 if (dx == 0.0) dx = 1.0; 2055 #if !defined(PETSC_USE_COMPLEX) 2056 if (dx < umin && dx >= 0.0) dx = umin; 2057 else if (dx < 0.0 && dx > -umin) dx = -umin; 2058 #else 2059 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2060 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2061 #endif 2062 dx *= epsilon; 2063 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2064 w3_array[col] += dx; 2065 } 2066 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2067 2068 /* 2069 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2070 */ 2071 2072 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2073 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2074 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2075 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 2076 2077 /* 2078 Loop over rows of vector, putting results into Jacobian matrix 2079 */ 2080 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2081 for (l=0; l<coloring->nrows[k]; l++) { 2082 row = coloring->rows[k][l]; 2083 col = coloring->columnsforrow[k][l]; 2084 y[row] *= vscale_array[vscaleforrow[k][l]]; 2085 srow = row + start; 2086 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2087 } 2088 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2089 } 2090 coloring->currentcolor = k; 2091 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2092 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2093 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2094 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2095 PetscFunctionReturn(0); 2096 } 2097 2098 #include "petscblaslapack.h" 2099 #undef __FUNCT__ 2100 #define __FUNCT__ "MatAXPY_SeqAIJ" 2101 PetscErrorCode MatAXPY_SeqAIJ(const PetscScalar a[],Mat X,Mat Y,MatStructure str) 2102 { 2103 PetscErrorCode ierr; 2104 PetscInt i; 2105 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2106 PetscBLASInt one=1,bnz = (PetscBLASInt)x->nz; 2107 2108 PetscFunctionBegin; 2109 if (str == SAME_NONZERO_PATTERN) { 2110 BLASaxpy_(&bnz,(PetscScalar*)a,x->a,&one,y->a,&one); 2111 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2112 if (y->xtoy && y->XtoY != X) { 2113 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2114 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2115 } 2116 if (!y->xtoy) { /* get xtoy */ 2117 ierr = MatAXPYGetxtoy_Private(X->m,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2118 y->XtoY = X; 2119 } 2120 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += (*a)*(x->a[i]); 2121 PetscLogInfo(0,"MatAXPY_SeqAIJ: ratio of nnz(X)/nnz(Y): %d/%d = %g\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz); 2122 } else { 2123 ierr = MatAXPY_Basic(a,X,Y,str);CHKERRQ(ierr); 2124 } 2125 PetscFunctionReturn(0); 2126 } 2127 2128 #undef __FUNCT__ 2129 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2130 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2131 { 2132 PetscFunctionBegin; 2133 PetscFunctionReturn(0); 2134 } 2135 2136 /* -------------------------------------------------------------------*/ 2137 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2138 MatGetRow_SeqAIJ, 2139 MatRestoreRow_SeqAIJ, 2140 MatMult_SeqAIJ, 2141 /* 4*/ MatMultAdd_SeqAIJ, 2142 MatMultTranspose_SeqAIJ, 2143 MatMultTransposeAdd_SeqAIJ, 2144 MatSolve_SeqAIJ, 2145 MatSolveAdd_SeqAIJ, 2146 MatSolveTranspose_SeqAIJ, 2147 /*10*/ MatSolveTransposeAdd_SeqAIJ, 2148 MatLUFactor_SeqAIJ, 2149 0, 2150 MatRelax_SeqAIJ, 2151 MatTranspose_SeqAIJ, 2152 /*15*/ MatGetInfo_SeqAIJ, 2153 MatEqual_SeqAIJ, 2154 MatGetDiagonal_SeqAIJ, 2155 MatDiagonalScale_SeqAIJ, 2156 MatNorm_SeqAIJ, 2157 /*20*/ 0, 2158 MatAssemblyEnd_SeqAIJ, 2159 MatCompress_SeqAIJ, 2160 MatSetOption_SeqAIJ, 2161 MatZeroEntries_SeqAIJ, 2162 /*25*/ MatZeroRows_SeqAIJ, 2163 MatLUFactorSymbolic_SeqAIJ, 2164 MatLUFactorNumeric_SeqAIJ, 2165 MatCholeskyFactorSymbolic_SeqAIJ, 2166 MatCholeskyFactorNumeric_SeqAIJ, 2167 /*30*/ MatSetUpPreallocation_SeqAIJ, 2168 MatILUFactorSymbolic_SeqAIJ, 2169 MatICCFactorSymbolic_SeqAIJ, 2170 MatGetArray_SeqAIJ, 2171 MatRestoreArray_SeqAIJ, 2172 /*35*/ MatDuplicate_SeqAIJ, 2173 0, 2174 0, 2175 MatILUFactor_SeqAIJ, 2176 0, 2177 /*40*/ MatAXPY_SeqAIJ, 2178 MatGetSubMatrices_SeqAIJ, 2179 MatIncreaseOverlap_SeqAIJ, 2180 MatGetValues_SeqAIJ, 2181 MatCopy_SeqAIJ, 2182 /*45*/ MatPrintHelp_SeqAIJ, 2183 MatScale_SeqAIJ, 2184 0, 2185 0, 2186 MatILUDTFactor_SeqAIJ, 2187 /*50*/ MatSetBlockSize_SeqAIJ, 2188 MatGetRowIJ_SeqAIJ, 2189 MatRestoreRowIJ_SeqAIJ, 2190 MatGetColumnIJ_SeqAIJ, 2191 MatRestoreColumnIJ_SeqAIJ, 2192 /*55*/ MatFDColoringCreate_SeqAIJ, 2193 0, 2194 0, 2195 MatPermute_SeqAIJ, 2196 0, 2197 /*60*/ 0, 2198 MatDestroy_SeqAIJ, 2199 MatView_SeqAIJ, 2200 MatGetPetscMaps_Petsc, 2201 0, 2202 /*65*/ 0, 2203 0, 2204 0, 2205 0, 2206 0, 2207 /*70*/ 0, 2208 0, 2209 MatSetColoring_SeqAIJ, 2210 #if defined(PETSC_HAVE_ADIC) 2211 MatSetValuesAdic_SeqAIJ, 2212 #else 2213 0, 2214 #endif 2215 MatSetValuesAdifor_SeqAIJ, 2216 /*75*/ MatFDColoringApply_SeqAIJ, 2217 0, 2218 0, 2219 0, 2220 0, 2221 /*80*/ 0, 2222 0, 2223 0, 2224 0, 2225 MatLoad_SeqAIJ, 2226 /*85*/ MatIsSymmetric_SeqAIJ, 2227 0, 2228 0, 2229 0, 2230 0, 2231 /*90*/ MatMatMult_SeqAIJ_SeqAIJ, 2232 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2233 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2234 MatPtAP_Basic, 2235 MatPtAPSymbolic_SeqAIJ, 2236 /*95*/ MatPtAPNumeric_SeqAIJ, 2237 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2238 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2239 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2240 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2241 /*100*/MatPtAPNumeric_SeqAIJ_SeqAIJ, 2242 0, 2243 0, 2244 }; 2245 2246 EXTERN_C_BEGIN 2247 #undef __FUNCT__ 2248 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2249 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2250 { 2251 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2252 PetscInt i,nz,n; 2253 2254 PetscFunctionBegin; 2255 2256 nz = aij->maxnz; 2257 n = mat->n; 2258 for (i=0; i<nz; i++) { 2259 aij->j[i] = indices[i]; 2260 } 2261 aij->nz = nz; 2262 for (i=0; i<n; i++) { 2263 aij->ilen[i] = aij->imax[i]; 2264 } 2265 2266 PetscFunctionReturn(0); 2267 } 2268 EXTERN_C_END 2269 2270 #undef __FUNCT__ 2271 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2272 /*@ 2273 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2274 in the matrix. 2275 2276 Input Parameters: 2277 + mat - the SeqAIJ matrix 2278 - indices - the column indices 2279 2280 Level: advanced 2281 2282 Notes: 2283 This can be called if you have precomputed the nonzero structure of the 2284 matrix and want to provide it to the matrix object to improve the performance 2285 of the MatSetValues() operation. 2286 2287 You MUST have set the correct numbers of nonzeros per row in the call to 2288 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2289 2290 MUST be called before any calls to MatSetValues(); 2291 2292 The indices should start with zero, not one. 2293 2294 @*/ 2295 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2296 { 2297 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2298 2299 PetscFunctionBegin; 2300 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2301 PetscValidPointer(indices,2); 2302 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2303 if (f) { 2304 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2305 } else { 2306 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2307 } 2308 PetscFunctionReturn(0); 2309 } 2310 2311 /* ----------------------------------------------------------------------------------------*/ 2312 2313 EXTERN_C_BEGIN 2314 #undef __FUNCT__ 2315 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2316 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 2317 { 2318 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2319 PetscErrorCode ierr; 2320 size_t nz = aij->i[mat->m]; 2321 2322 PetscFunctionBegin; 2323 if (aij->nonew != 1) { 2324 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2325 } 2326 2327 /* allocate space for values if not already there */ 2328 if (!aij->saved_values) { 2329 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2330 } 2331 2332 /* copy values over */ 2333 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2334 PetscFunctionReturn(0); 2335 } 2336 EXTERN_C_END 2337 2338 #undef __FUNCT__ 2339 #define __FUNCT__ "MatStoreValues" 2340 /*@ 2341 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2342 example, reuse of the linear part of a Jacobian, while recomputing the 2343 nonlinear portion. 2344 2345 Collect on Mat 2346 2347 Input Parameters: 2348 . mat - the matrix (currently only AIJ matrices support this option) 2349 2350 Level: advanced 2351 2352 Common Usage, with SNESSolve(): 2353 $ Create Jacobian matrix 2354 $ Set linear terms into matrix 2355 $ Apply boundary conditions to matrix, at this time matrix must have 2356 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2357 $ boundary conditions again will not change the nonzero structure 2358 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2359 $ ierr = MatStoreValues(mat); 2360 $ Call SNESSetJacobian() with matrix 2361 $ In your Jacobian routine 2362 $ ierr = MatRetrieveValues(mat); 2363 $ Set nonlinear terms in matrix 2364 2365 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2366 $ // build linear portion of Jacobian 2367 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2368 $ ierr = MatStoreValues(mat); 2369 $ loop over nonlinear iterations 2370 $ ierr = MatRetrieveValues(mat); 2371 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2372 $ // call MatAssemblyBegin/End() on matrix 2373 $ Solve linear system with Jacobian 2374 $ endloop 2375 2376 Notes: 2377 Matrix must already be assemblied before calling this routine 2378 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2379 calling this routine. 2380 2381 .seealso: MatRetrieveValues() 2382 2383 @*/ 2384 PetscErrorCode MatStoreValues(Mat mat) 2385 { 2386 PetscErrorCode ierr,(*f)(Mat); 2387 2388 PetscFunctionBegin; 2389 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2390 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2391 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2392 2393 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2394 if (f) { 2395 ierr = (*f)(mat);CHKERRQ(ierr); 2396 } else { 2397 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2398 } 2399 PetscFunctionReturn(0); 2400 } 2401 2402 EXTERN_C_BEGIN 2403 #undef __FUNCT__ 2404 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2405 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 2406 { 2407 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2408 PetscErrorCode ierr; 2409 PetscInt nz = aij->i[mat->m]; 2410 2411 PetscFunctionBegin; 2412 if (aij->nonew != 1) { 2413 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2414 } 2415 if (!aij->saved_values) { 2416 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2417 } 2418 /* copy values over */ 2419 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2420 PetscFunctionReturn(0); 2421 } 2422 EXTERN_C_END 2423 2424 #undef __FUNCT__ 2425 #define __FUNCT__ "MatRetrieveValues" 2426 /*@ 2427 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2428 example, reuse of the linear part of a Jacobian, while recomputing the 2429 nonlinear portion. 2430 2431 Collect on Mat 2432 2433 Input Parameters: 2434 . mat - the matrix (currently on AIJ matrices support this option) 2435 2436 Level: advanced 2437 2438 .seealso: MatStoreValues() 2439 2440 @*/ 2441 PetscErrorCode MatRetrieveValues(Mat mat) 2442 { 2443 PetscErrorCode ierr,(*f)(Mat); 2444 2445 PetscFunctionBegin; 2446 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2447 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2448 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2449 2450 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2451 if (f) { 2452 ierr = (*f)(mat);CHKERRQ(ierr); 2453 } else { 2454 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2455 } 2456 PetscFunctionReturn(0); 2457 } 2458 2459 2460 /* --------------------------------------------------------------------------------*/ 2461 #undef __FUNCT__ 2462 #define __FUNCT__ "MatCreateSeqAIJ" 2463 /*@C 2464 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2465 (the default parallel PETSc format). For good matrix assembly performance 2466 the user should preallocate the matrix storage by setting the parameter nz 2467 (or the array nnz). By setting these parameters accurately, performance 2468 during matrix assembly can be increased by more than a factor of 50. 2469 2470 Collective on MPI_Comm 2471 2472 Input Parameters: 2473 + comm - MPI communicator, set to PETSC_COMM_SELF 2474 . m - number of rows 2475 . n - number of columns 2476 . nz - number of nonzeros per row (same for all rows) 2477 - nnz - array containing the number of nonzeros in the various rows 2478 (possibly different for each row) or PETSC_NULL 2479 2480 Output Parameter: 2481 . A - the matrix 2482 2483 Notes: 2484 If nnz is given then nz is ignored 2485 2486 The AIJ format (also called the Yale sparse matrix format or 2487 compressed row storage), is fully compatible with standard Fortran 77 2488 storage. That is, the stored row and column indices can begin at 2489 either one (as in Fortran) or zero. See the users' manual for details. 2490 2491 Specify the preallocated storage with either nz or nnz (not both). 2492 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2493 allocation. For large problems you MUST preallocate memory or you 2494 will get TERRIBLE performance, see the users' manual chapter on matrices. 2495 2496 By default, this format uses inodes (identical nodes) when possible, to 2497 improve numerical efficiency of matrix-vector products and solves. We 2498 search for consecutive rows with the same nonzero structure, thereby 2499 reusing matrix information to achieve increased efficiency. 2500 2501 Options Database Keys: 2502 + -mat_no_inode - Do not use inodes 2503 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2504 - -mat_aij_oneindex - Internally use indexing starting at 1 2505 rather than 0. Note that when calling MatSetValues(), 2506 the user still MUST index entries starting at 0! 2507 2508 Level: intermediate 2509 2510 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2511 2512 @*/ 2513 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2514 { 2515 PetscErrorCode ierr; 2516 2517 PetscFunctionBegin; 2518 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2519 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2520 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2521 PetscFunctionReturn(0); 2522 } 2523 2524 #define SKIP_ALLOCATION -4 2525 2526 #undef __FUNCT__ 2527 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2528 /*@C 2529 MatSeqAIJSetPreallocation - For good matrix assembly performance 2530 the user should preallocate the matrix storage by setting the parameter nz 2531 (or the array nnz). By setting these parameters accurately, performance 2532 during matrix assembly can be increased by more than a factor of 50. 2533 2534 Collective on MPI_Comm 2535 2536 Input Parameters: 2537 + comm - MPI communicator, set to PETSC_COMM_SELF 2538 . m - number of rows 2539 . n - number of columns 2540 . nz - number of nonzeros per row (same for all rows) 2541 - nnz - array containing the number of nonzeros in the various rows 2542 (possibly different for each row) or PETSC_NULL 2543 2544 Output Parameter: 2545 . A - the matrix 2546 2547 Notes: 2548 If nnz is given then nz is ignored 2549 2550 The AIJ format (also called the Yale sparse matrix format or 2551 compressed row storage), is fully compatible with standard Fortran 77 2552 storage. That is, the stored row and column indices can begin at 2553 either one (as in Fortran) or zero. See the users' manual for details. 2554 2555 Specify the preallocated storage with either nz or nnz (not both). 2556 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2557 allocation. For large problems you MUST preallocate memory or you 2558 will get TERRIBLE performance, see the users' manual chapter on matrices. 2559 2560 By default, this format uses inodes (identical nodes) when possible, to 2561 improve numerical efficiency of matrix-vector products and solves. We 2562 search for consecutive rows with the same nonzero structure, thereby 2563 reusing matrix information to achieve increased efficiency. 2564 2565 Options Database Keys: 2566 + -mat_no_inode - Do not use inodes 2567 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2568 - -mat_aij_oneindex - Internally use indexing starting at 1 2569 rather than 0. Note that when calling MatSetValues(), 2570 the user still MUST index entries starting at 0! 2571 2572 Level: intermediate 2573 2574 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2575 2576 @*/ 2577 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2578 { 2579 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2580 2581 PetscFunctionBegin; 2582 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2583 if (f) { 2584 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2585 } 2586 PetscFunctionReturn(0); 2587 } 2588 2589 EXTERN_C_BEGIN 2590 #undef __FUNCT__ 2591 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2592 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2593 { 2594 Mat_SeqAIJ *b; 2595 size_t len = 0; 2596 PetscTruth skipallocation = PETSC_FALSE; 2597 PetscErrorCode ierr; 2598 PetscInt i; 2599 2600 PetscFunctionBegin; 2601 2602 if (nz == SKIP_ALLOCATION) { 2603 skipallocation = PETSC_TRUE; 2604 nz = 0; 2605 } 2606 2607 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2608 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2609 if (nnz) { 2610 for (i=0; i<B->m; i++) { 2611 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2612 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); 2613 } 2614 } 2615 2616 B->preallocated = PETSC_TRUE; 2617 b = (Mat_SeqAIJ*)B->data; 2618 2619 ierr = PetscMalloc((B->m+1)*sizeof(PetscInt),&b->imax);CHKERRQ(ierr); 2620 if (!nnz) { 2621 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2622 else if (nz <= 0) nz = 1; 2623 for (i=0; i<B->m; i++) b->imax[i] = nz; 2624 nz = nz*B->m; 2625 } else { 2626 nz = 0; 2627 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2628 } 2629 2630 if (!skipallocation) { 2631 /* allocate the matrix space */ 2632 len = ((size_t) nz)*(sizeof(PetscInt) + sizeof(PetscScalar)) + (B->m+1)*sizeof(PetscInt); 2633 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2634 b->j = (PetscInt*)(b->a + nz); 2635 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2636 b->i = b->j + nz; 2637 b->i[0] = 0; 2638 for (i=1; i<B->m+1; i++) { 2639 b->i[i] = b->i[i-1] + b->imax[i-1]; 2640 } 2641 b->singlemalloc = PETSC_TRUE; 2642 b->freedata = PETSC_TRUE; 2643 } else { 2644 b->freedata = PETSC_FALSE; 2645 } 2646 2647 /* b->ilen will count nonzeros in each row so far. */ 2648 ierr = PetscMalloc((B->m+1)*sizeof(PetscInt),&b->ilen);CHKERRQ(ierr); 2649 ierr = PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(PetscInt)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2650 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2651 2652 b->nz = 0; 2653 b->maxnz = nz; 2654 B->info.nz_unneeded = (double)b->maxnz; 2655 PetscFunctionReturn(0); 2656 } 2657 EXTERN_C_END 2658 2659 /*MC 2660 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 2661 based on compressed sparse row format. 2662 2663 Options Database Keys: 2664 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 2665 2666 Level: beginner 2667 2668 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 2669 M*/ 2670 2671 EXTERN_C_BEGIN 2672 #undef __FUNCT__ 2673 #define __FUNCT__ "MatCreate_SeqAIJ" 2674 PetscErrorCode MatCreate_SeqAIJ(Mat B) 2675 { 2676 Mat_SeqAIJ *b; 2677 PetscErrorCode ierr; 2678 PetscMPIInt size; 2679 2680 PetscFunctionBegin; 2681 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2682 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2683 2684 B->m = B->M = PetscMax(B->m,B->M); 2685 B->n = B->N = PetscMax(B->n,B->N); 2686 2687 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2688 B->data = (void*)b; 2689 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2690 B->factor = 0; 2691 B->lupivotthreshold = 1.0; 2692 B->mapping = 0; 2693 ierr = PetscOptionsGetReal(B->prefix,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2694 ierr = PetscOptionsHasName(B->prefix,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2695 b->row = 0; 2696 b->col = 0; 2697 b->icol = 0; 2698 b->reallocs = 0; 2699 b->lu_shift = PETSC_FALSE; 2700 2701 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2702 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2703 2704 b->sorted = PETSC_FALSE; 2705 b->ignorezeroentries = PETSC_FALSE; 2706 b->roworiented = PETSC_TRUE; 2707 b->nonew = 0; 2708 b->diag = 0; 2709 b->solve_work = 0; 2710 B->spptr = 0; 2711 b->saved_values = 0; 2712 b->idiag = 0; 2713 b->ssor = 0; 2714 b->keepzeroedrows = PETSC_FALSE; 2715 b->xtoy = 0; 2716 b->XtoY = 0; 2717 b->compressedrow.use = PETSC_FALSE; 2718 b->compressedrow.nrows = B->m; 2719 b->compressedrow.i = PETSC_NULL; 2720 b->compressedrow.rindex = PETSC_NULL; 2721 b->compressedrow.checked = PETSC_FALSE; 2722 B->same_nonzero = PETSC_FALSE; 2723 2724 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2725 2726 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2727 "MatSeqAIJSetColumnIndices_SeqAIJ", 2728 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2729 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2730 "MatStoreValues_SeqAIJ", 2731 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2732 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2733 "MatRetrieveValues_SeqAIJ", 2734 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2735 #if !defined(PETSC_USE_64BIT_INT) 2736 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2737 "MatConvert_SeqAIJ_SeqSBAIJ", 2738 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2739 #endif 2740 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2741 "MatConvert_SeqAIJ_SeqBAIJ", 2742 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2743 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 2744 "MatIsTranspose_SeqAIJ", 2745 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 2746 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2747 "MatSeqAIJSetPreallocation_SeqAIJ", 2748 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2749 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 2750 "MatReorderForNonzeroDiagonal_SeqAIJ", 2751 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 2752 ierr = MatCreate_Inode(B);CHKERRQ(ierr); 2753 PetscFunctionReturn(0); 2754 } 2755 EXTERN_C_END 2756 2757 #undef __FUNCT__ 2758 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2759 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2760 { 2761 Mat C; 2762 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2763 PetscErrorCode ierr; 2764 PetscInt i,m = A->m; 2765 size_t len; 2766 2767 PetscFunctionBegin; 2768 *B = 0; 2769 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2770 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 2771 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2772 2773 c = (Mat_SeqAIJ*)C->data; 2774 2775 C->factor = A->factor; 2776 C->lupivotthreshold = A->lupivotthreshold; 2777 2778 c->row = 0; 2779 c->col = 0; 2780 c->icol = 0; 2781 c->reallocs = 0; 2782 c->lu_shift = PETSC_FALSE; 2783 2784 C->assembled = PETSC_TRUE; 2785 2786 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->imax);CHKERRQ(ierr); 2787 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->ilen);CHKERRQ(ierr); 2788 for (i=0; i<m; i++) { 2789 c->imax[i] = a->imax[i]; 2790 c->ilen[i] = a->ilen[i]; 2791 } 2792 2793 /* allocate the matrix space */ 2794 c->singlemalloc = PETSC_TRUE; 2795 len = ((size_t) (m+1))*sizeof(PetscInt)+(a->i[m])*(sizeof(PetscScalar)+sizeof(PetscInt)); 2796 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2797 c->j = (PetscInt*)(c->a + a->i[m] ); 2798 c->i = c->j + a->i[m]; 2799 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2800 if (m > 0) { 2801 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 2802 if (cpvalues == MAT_COPY_VALUES) { 2803 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2804 } else { 2805 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2806 } 2807 } 2808 2809 ierr = PetscLogObjectMemory(C,len+2*(m+1)*sizeof(PetscInt)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2810 c->sorted = a->sorted; 2811 c->ignorezeroentries = a->ignorezeroentries; 2812 c->roworiented = a->roworiented; 2813 c->nonew = a->nonew; 2814 if (a->diag) { 2815 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2816 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2817 for (i=0; i<m; i++) { 2818 c->diag[i] = a->diag[i]; 2819 } 2820 } else c->diag = 0; 2821 c->solve_work = 0; 2822 c->saved_values = 0; 2823 c->idiag = 0; 2824 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2825 c->ssor = 0; 2826 c->keepzeroedrows = a->keepzeroedrows; 2827 c->freedata = PETSC_TRUE; 2828 c->xtoy = 0; 2829 c->XtoY = 0; 2830 2831 c->nz = a->nz; 2832 c->maxnz = a->maxnz; 2833 C->preallocated = PETSC_TRUE; 2834 2835 c->compressedrow.use = a->compressedrow.use; 2836 c->compressedrow.nrows = a->compressedrow.nrows; 2837 c->compressedrow.checked = a->compressedrow.checked; 2838 if ( a->compressedrow.checked && a->compressedrow.use){ 2839 i = a->compressedrow.nrows; 2840 ierr = PetscMalloc((2*i+1)*sizeof(PetscInt),&c->compressedrow.i);CHKERRQ(ierr); 2841 c->compressedrow.rindex = c->compressedrow.i + i + 1; 2842 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 2843 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 2844 } else { 2845 c->compressedrow.use = PETSC_FALSE; 2846 c->compressedrow.i = PETSC_NULL; 2847 c->compressedrow.rindex = PETSC_NULL; 2848 } 2849 C->same_nonzero = A->same_nonzero; 2850 2851 ierr = MatDuplicate_Inode(A,cpvalues,&C);CHKERRQ(ierr); 2852 2853 *B = C; 2854 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2855 PetscFunctionReturn(0); 2856 } 2857 2858 #undef __FUNCT__ 2859 #define __FUNCT__ "MatLoad_SeqAIJ" 2860 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer,const MatType type,Mat *A) 2861 { 2862 Mat_SeqAIJ *a; 2863 Mat B; 2864 PetscErrorCode ierr; 2865 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 2866 int fd; 2867 PetscMPIInt size; 2868 MPI_Comm comm; 2869 2870 PetscFunctionBegin; 2871 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2872 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2873 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2874 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2875 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2876 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2877 M = header[1]; N = header[2]; nz = header[3]; 2878 2879 if (nz < 0) { 2880 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2881 } 2882 2883 /* read in row lengths */ 2884 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2885 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2886 2887 /* check if sum of rowlengths is same as nz */ 2888 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 2889 if (sum != nz) SETERRQ2(PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 2890 2891 /* create our matrix */ 2892 ierr = MatCreate(comm,PETSC_DECIDE,PETSC_DECIDE,M,N,&B);CHKERRQ(ierr); 2893 ierr = MatSetType(B,type);CHKERRQ(ierr); 2894 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 2895 a = (Mat_SeqAIJ*)B->data; 2896 2897 /* read in column indices and adjust for Fortran indexing*/ 2898 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2899 2900 /* read in nonzero values */ 2901 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2902 2903 /* set matrix "i" values */ 2904 a->i[0] = 0; 2905 for (i=1; i<= M; i++) { 2906 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2907 a->ilen[i-1] = rowlengths[i-1]; 2908 } 2909 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2910 2911 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2912 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2913 *A = B; 2914 PetscFunctionReturn(0); 2915 } 2916 2917 #undef __FUNCT__ 2918 #define __FUNCT__ "MatEqual_SeqAIJ" 2919 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2920 { 2921 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2922 PetscErrorCode ierr; 2923 2924 PetscFunctionBegin; 2925 /* If the matrix dimensions are not equal,or no of nonzeros */ 2926 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)) { 2927 *flg = PETSC_FALSE; 2928 PetscFunctionReturn(0); 2929 } 2930 2931 /* if the a->i are the same */ 2932 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 2933 if (!*flg) PetscFunctionReturn(0); 2934 2935 /* if a->j are the same */ 2936 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 2937 if (!*flg) PetscFunctionReturn(0); 2938 2939 /* if a->a are the same */ 2940 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2941 2942 PetscFunctionReturn(0); 2943 2944 } 2945 2946 #undef __FUNCT__ 2947 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2948 /*@C 2949 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2950 provided by the user. 2951 2952 Coolective on MPI_Comm 2953 2954 Input Parameters: 2955 + comm - must be an MPI communicator of size 1 2956 . m - number of rows 2957 . n - number of columns 2958 . i - row indices 2959 . j - column indices 2960 - a - matrix values 2961 2962 Output Parameter: 2963 . mat - the matrix 2964 2965 Level: intermediate 2966 2967 Notes: 2968 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2969 once the matrix is destroyed 2970 2971 You cannot set new nonzero locations into this matrix, that will generate an error. 2972 2973 The i and j indices are 0 based 2974 2975 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2976 2977 @*/ 2978 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2979 { 2980 PetscErrorCode ierr; 2981 PetscInt ii; 2982 Mat_SeqAIJ *aij; 2983 2984 PetscFunctionBegin; 2985 ierr = MatCreate(comm,m,n,m,n,mat);CHKERRQ(ierr); 2986 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 2987 ierr = MatSeqAIJSetPreallocation(*mat,SKIP_ALLOCATION,0);CHKERRQ(ierr); 2988 aij = (Mat_SeqAIJ*)(*mat)->data; 2989 2990 if (i[0] != 0) { 2991 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2992 } 2993 aij->i = i; 2994 aij->j = j; 2995 aij->a = a; 2996 aij->singlemalloc = PETSC_FALSE; 2997 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2998 aij->freedata = PETSC_FALSE; 2999 3000 for (ii=0; ii<m; ii++) { 3001 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3002 #if defined(PETSC_USE_DEBUG) 3003 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]); 3004 #endif 3005 } 3006 #if defined(PETSC_USE_DEBUG) 3007 for (ii=0; ii<aij->i[m]; ii++) { 3008 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3009 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3010 } 3011 #endif 3012 3013 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3014 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3015 PetscFunctionReturn(0); 3016 } 3017 3018 #undef __FUNCT__ 3019 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3020 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3021 { 3022 PetscErrorCode ierr; 3023 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3024 3025 PetscFunctionBegin; 3026 if (coloring->ctype == IS_COLORING_LOCAL) { 3027 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3028 a->coloring = coloring; 3029 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3030 PetscInt i,*larray; 3031 ISColoring ocoloring; 3032 ISColoringValue *colors; 3033 3034 /* set coloring for diagonal portion */ 3035 ierr = PetscMalloc((A->n+1)*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3036 for (i=0; i<A->n; i++) { 3037 larray[i] = i; 3038 } 3039 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3040 ierr = PetscMalloc((A->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3041 for (i=0; i<A->n; i++) { 3042 colors[i] = coloring->colors[larray[i]]; 3043 } 3044 ierr = PetscFree(larray);CHKERRQ(ierr); 3045 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 3046 a->coloring = ocoloring; 3047 } 3048 PetscFunctionReturn(0); 3049 } 3050 3051 #if defined(PETSC_HAVE_ADIC) 3052 EXTERN_C_BEGIN 3053 #include "adic/ad_utils.h" 3054 EXTERN_C_END 3055 3056 #undef __FUNCT__ 3057 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3058 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3059 { 3060 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3061 PetscInt m = A->m,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3062 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3063 ISColoringValue *color; 3064 3065 PetscFunctionBegin; 3066 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3067 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3068 color = a->coloring->colors; 3069 /* loop over rows */ 3070 for (i=0; i<m; i++) { 3071 nz = ii[i+1] - ii[i]; 3072 /* loop over columns putting computed value into matrix */ 3073 for (j=0; j<nz; j++) { 3074 *v++ = values[color[*jj++]]; 3075 } 3076 values += nlen; /* jump to next row of derivatives */ 3077 } 3078 PetscFunctionReturn(0); 3079 } 3080 #endif 3081 3082 #undef __FUNCT__ 3083 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3084 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3085 { 3086 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3087 PetscInt m = A->m,*ii = a->i,*jj = a->j,nz,i,j; 3088 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3089 ISColoringValue *color; 3090 3091 PetscFunctionBegin; 3092 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3093 color = a->coloring->colors; 3094 /* loop over rows */ 3095 for (i=0; i<m; i++) { 3096 nz = ii[i+1] - ii[i]; 3097 /* loop over columns putting computed value into matrix */ 3098 for (j=0; j<nz; j++) { 3099 *v++ = values[color[*jj++]]; 3100 } 3101 values += nl; /* jump to next row of derivatives */ 3102 } 3103 PetscFunctionReturn(0); 3104 } 3105 3106