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