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