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