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