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 PetscLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 807 break; 808 case MAT_NO_NEW_DIAGONALS: 809 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 810 case MAT_INODE_LIMIT_1: 811 a->inode.limit = 1; 812 break; 813 case MAT_INODE_LIMIT_2: 814 a->inode.limit = 2; 815 break; 816 case MAT_INODE_LIMIT_3: 817 a->inode.limit = 3; 818 break; 819 case MAT_INODE_LIMIT_4: 820 a->inode.limit = 4; 821 break; 822 case MAT_INODE_LIMIT_5: 823 a->inode.limit = 5; 824 break; 825 default: 826 SETERRQ(PETSC_ERR_SUP,"unknown option"); 827 } 828 PetscFunctionReturn(0); 829 } 830 831 #undef __FUNCT__ 832 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 833 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 834 { 835 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 836 int i,j,n,shift = a->indexshift,ierr; 837 PetscScalar *x,zero = 0.0; 838 839 PetscFunctionBegin; 840 ierr = VecSet(&zero,v);CHKERRQ(ierr); 841 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 842 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 843 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 844 for (i=0; i<A->m; i++) { 845 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 846 if (a->j[j]+shift == i) { 847 x[i] = a->a[j]; 848 break; 849 } 850 } 851 } 852 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 853 PetscFunctionReturn(0); 854 } 855 856 857 #undef __FUNCT__ 858 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 859 int MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 860 { 861 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 862 PetscScalar *x,*y; 863 int ierr,m = A->m,shift = a->indexshift; 864 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 865 PetscScalar *v,alpha; 866 int n,i,*idx; 867 #endif 868 869 PetscFunctionBegin; 870 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 871 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 872 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 873 y = y + shift; /* shift for Fortran start by 1 indexing */ 874 875 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 876 fortranmulttransposeaddaij_(&m,x,a->i,a->j+shift,a->a+shift,y); 877 #else 878 for (i=0; i<m; i++) { 879 idx = a->j + a->i[i] + shift; 880 v = a->a + a->i[i] + shift; 881 n = a->i[i+1] - a->i[i]; 882 alpha = x[i]; 883 while (n-->0) {y[*idx++] += alpha * *v++;} 884 } 885 #endif 886 PetscLogFlops(2*a->nz); 887 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 888 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 889 PetscFunctionReturn(0); 890 } 891 892 #undef __FUNCT__ 893 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 894 int MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 895 { 896 PetscScalar zero = 0.0; 897 int ierr; 898 899 PetscFunctionBegin; 900 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 901 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 902 PetscFunctionReturn(0); 903 } 904 905 906 #undef __FUNCT__ 907 #define __FUNCT__ "MatMult_SeqAIJ" 908 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 909 { 910 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 911 PetscScalar *x,*y,*v; 912 int ierr,m = A->m,*idx,shift = a->indexshift,*ii; 913 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 914 int n,i,jrow,j; 915 PetscScalar sum; 916 #endif 917 918 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 919 #pragma disjoint(*x,*y,*v) 920 #endif 921 922 PetscFunctionBegin; 923 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 924 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 925 x = x + shift; /* shift for Fortran start by 1 indexing */ 926 idx = a->j; 927 v = a->a; 928 ii = a->i; 929 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 930 fortranmultaij_(&m,x,ii,idx+shift,v+shift,y); 931 #else 932 v += shift; /* shift for Fortran start by 1 indexing */ 933 idx += shift; 934 for (i=0; i<m; i++) { 935 jrow = ii[i]; 936 n = ii[i+1] - jrow; 937 sum = 0.0; 938 for (j=0; j<n; j++) { 939 sum += v[jrow]*x[idx[jrow]]; jrow++; 940 } 941 y[i] = sum; 942 } 943 #endif 944 PetscLogFlops(2*a->nz - m); 945 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 946 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 947 PetscFunctionReturn(0); 948 } 949 950 #undef __FUNCT__ 951 #define __FUNCT__ "MatMultAdd_SeqAIJ" 952 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 953 { 954 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 955 PetscScalar *x,*y,*z,*v; 956 int ierr,m = A->m,*idx,shift = a->indexshift,*ii; 957 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 958 int n,i,jrow,j; 959 PetscScalar sum; 960 #endif 961 962 PetscFunctionBegin; 963 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 964 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 965 if (zz != yy) { 966 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 967 } else { 968 z = y; 969 } 970 x = x + shift; /* shift for Fortran start by 1 indexing */ 971 idx = a->j; 972 v = a->a; 973 ii = a->i; 974 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 975 fortranmultaddaij_(&m,x,ii,idx+shift,v+shift,y,z); 976 #else 977 v += shift; /* shift for Fortran start by 1 indexing */ 978 idx += shift; 979 for (i=0; i<m; i++) { 980 jrow = ii[i]; 981 n = ii[i+1] - jrow; 982 sum = y[i]; 983 for (j=0; j<n; j++) { 984 sum += v[jrow]*x[idx[jrow]]; jrow++; 985 } 986 z[i] = sum; 987 } 988 #endif 989 PetscLogFlops(2*a->nz); 990 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 991 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 992 if (zz != yy) { 993 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 994 } 995 PetscFunctionReturn(0); 996 } 997 998 /* 999 Adds diagonal pointers to sparse matrix structure. 1000 */ 1001 #undef __FUNCT__ 1002 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1003 int MatMarkDiagonal_SeqAIJ(Mat A) 1004 { 1005 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1006 int i,j,*diag,m = A->m,shift = a->indexshift,ierr; 1007 1008 PetscFunctionBegin; 1009 if (a->diag) PetscFunctionReturn(0); 1010 1011 ierr = PetscMalloc((m+1)*sizeof(int),&diag);CHKERRQ(ierr); 1012 PetscLogObjectMemory(A,(m+1)*sizeof(int)); 1013 for (i=0; i<A->m; i++) { 1014 diag[i] = a->i[i+1]; 1015 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 1016 if (a->j[j]+shift == i) { 1017 diag[i] = j - shift; 1018 break; 1019 } 1020 } 1021 } 1022 a->diag = diag; 1023 PetscFunctionReturn(0); 1024 } 1025 1026 /* 1027 Checks for missing diagonals 1028 */ 1029 #undef __FUNCT__ 1030 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1031 int MatMissingDiagonal_SeqAIJ(Mat A) 1032 { 1033 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1034 int *diag,*jj = a->j,i,shift = a->indexshift,ierr; 1035 1036 PetscFunctionBegin; 1037 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1038 diag = a->diag; 1039 for (i=0; i<A->m; i++) { 1040 if (jj[diag[i]+shift] != i-shift) { 1041 SETERRQ1(1,"Matrix is missing diagonal number %d",i); 1042 } 1043 } 1044 PetscFunctionReturn(0); 1045 } 1046 1047 #undef __FUNCT__ 1048 #define __FUNCT__ "MatRelax_SeqAIJ" 1049 int MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,int its,int lits,Vec xx) 1050 { 1051 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1052 PetscScalar *x,*b,*bs, d,*xs,sum,*v = a->a,*t=0,scale,*ts,*xb,*idiag=0; 1053 int ierr,*idx,*diag,n = A->n,m = A->m,i,shift = a->indexshift; 1054 1055 PetscFunctionBegin; 1056 its = its*lits; 1057 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %d and local its %d both positive",its,lits); 1058 1059 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1060 if (xx != bb) { 1061 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 1062 } else { 1063 b = x; 1064 } 1065 1066 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 1067 diag = a->diag; 1068 xs = x + shift; /* shifted by one for index start of a or a->j*/ 1069 if (flag == SOR_APPLY_UPPER) { 1070 /* apply (U + D/omega) to the vector */ 1071 bs = b + shift; 1072 for (i=0; i<m; i++) { 1073 d = fshift + a->a[diag[i] + shift]; 1074 n = a->i[i+1] - diag[i] - 1; 1075 PetscLogFlops(2*n-1); 1076 idx = a->j + diag[i] + (!shift); 1077 v = a->a + diag[i] + (!shift); 1078 sum = b[i]*d/omega; 1079 SPARSEDENSEDOT(sum,bs,v,idx,n); 1080 x[i] = sum; 1081 } 1082 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1083 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1084 PetscFunctionReturn(0); 1085 } 1086 1087 /* setup workspace for Eisenstat */ 1088 if (flag & SOR_EISENSTAT) { 1089 if (!a->idiag) { 1090 ierr = PetscMalloc(2*m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 1091 a->ssor = a->idiag + m; 1092 v = a->a; 1093 for (i=0; i<m; i++) { a->idiag[i] = 1.0/v[diag[i]];} 1094 } 1095 t = a->ssor; 1096 idiag = a->idiag; 1097 } 1098 /* Let A = L + U + D; where L is lower trianglar, 1099 U is upper triangular, E is diagonal; This routine applies 1100 1101 (L + E)^{-1} A (U + E)^{-1} 1102 1103 to a vector efficiently using Eisenstat's trick. This is for 1104 the case of SSOR preconditioner, so E is D/omega where omega 1105 is the relaxation factor. 1106 */ 1107 1108 if (flag == SOR_APPLY_LOWER) { 1109 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1110 } else if ((flag & SOR_EISENSTAT) && omega == 1.0 && shift == 0 && fshift == 0.0) { 1111 /* special case for omega = 1.0 saves flops and some integer ops */ 1112 PetscScalar *v2; 1113 1114 v2 = a->a; 1115 /* x = (E + U)^{-1} b */ 1116 for (i=m-1; i>=0; i--) { 1117 n = a->i[i+1] - diag[i] - 1; 1118 idx = a->j + diag[i] + 1; 1119 v = a->a + diag[i] + 1; 1120 sum = b[i]; 1121 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1122 x[i] = sum*idiag[i]; 1123 1124 /* t = b - (2*E - D)x */ 1125 t[i] = b[i] - (v2[diag[i]])*x[i]; 1126 } 1127 1128 /* t = (E + L)^{-1}t */ 1129 diag = a->diag; 1130 for (i=0; i<m; i++) { 1131 n = diag[i] - a->i[i]; 1132 idx = a->j + a->i[i]; 1133 v = a->a + a->i[i]; 1134 sum = t[i]; 1135 SPARSEDENSEMDOT(sum,t,v,idx,n); 1136 t[i] = sum*idiag[i]; 1137 1138 /* x = x + t */ 1139 x[i] += t[i]; 1140 } 1141 1142 PetscLogFlops(3*m-1 + 2*a->nz); 1143 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1144 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1145 PetscFunctionReturn(0); 1146 } else if (flag & SOR_EISENSTAT) { 1147 /* Let A = L + U + D; where L is lower trianglar, 1148 U is upper triangular, E is diagonal; This routine applies 1149 1150 (L + E)^{-1} A (U + E)^{-1} 1151 1152 to a vector efficiently using Eisenstat's trick. This is for 1153 the case of SSOR preconditioner, so E is D/omega where omega 1154 is the relaxation factor. 1155 */ 1156 scale = (2.0/omega) - 1.0; 1157 1158 /* x = (E + U)^{-1} b */ 1159 for (i=m-1; i>=0; i--) { 1160 d = fshift + a->a[diag[i] + shift]; 1161 n = a->i[i+1] - diag[i] - 1; 1162 idx = a->j + diag[i] + (!shift); 1163 v = a->a + diag[i] + (!shift); 1164 sum = b[i]; 1165 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1166 x[i] = omega*(sum/d); 1167 } 1168 1169 /* t = b - (2*E - D)x */ 1170 v = a->a; 1171 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 1172 1173 /* t = (E + L)^{-1}t */ 1174 ts = t + shift; /* shifted by one for index start of a or a->j*/ 1175 diag = a->diag; 1176 for (i=0; i<m; i++) { 1177 d = fshift + a->a[diag[i]+shift]; 1178 n = diag[i] - a->i[i]; 1179 idx = a->j + a->i[i] + shift; 1180 v = a->a + a->i[i] + shift; 1181 sum = t[i]; 1182 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1183 t[i] = omega*(sum/d); 1184 /* x = x + t */ 1185 x[i] += t[i]; 1186 } 1187 1188 PetscLogFlops(6*m-1 + 2*a->nz); 1189 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1190 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1191 PetscFunctionReturn(0); 1192 } 1193 if (flag & SOR_ZERO_INITIAL_GUESS) { 1194 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1195 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1196 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,diag,a->a,b); 1197 #else 1198 for (i=0; i<m; i++) { 1199 d = fshift + a->a[diag[i]+shift]; 1200 n = diag[i] - a->i[i]; 1201 PetscLogFlops(2*n-1); 1202 idx = a->j + a->i[i] + shift; 1203 v = a->a + a->i[i] + shift; 1204 sum = b[i]; 1205 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1206 x[i] = omega*(sum/d); 1207 } 1208 #endif 1209 xb = x; 1210 } else xb = b; 1211 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1212 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1213 for (i=0; i<m; i++) { 1214 x[i] *= a->a[diag[i]+shift]; 1215 } 1216 PetscLogFlops(m); 1217 } 1218 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1219 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1220 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,diag,a->a,xb); 1221 #else 1222 for (i=m-1; i>=0; i--) { 1223 d = fshift + a->a[diag[i] + shift]; 1224 n = a->i[i+1] - diag[i] - 1; 1225 PetscLogFlops(2*n-1); 1226 idx = a->j + diag[i] + (!shift); 1227 v = a->a + diag[i] + (!shift); 1228 sum = xb[i]; 1229 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1230 x[i] = omega*(sum/d); 1231 } 1232 #endif 1233 } 1234 its--; 1235 } 1236 while (its--) { 1237 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1238 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1239 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,diag,a->a,b); 1240 #else 1241 for (i=0; i<m; i++) { 1242 d = fshift + a->a[diag[i]+shift]; 1243 n = a->i[i+1] - a->i[i]; 1244 PetscLogFlops(2*n-1); 1245 idx = a->j + a->i[i] + shift; 1246 v = a->a + a->i[i] + shift; 1247 sum = b[i]; 1248 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1249 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1250 } 1251 #endif 1252 } 1253 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1254 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1255 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,diag,a->a,b); 1256 #else 1257 for (i=m-1; i>=0; i--) { 1258 d = fshift + a->a[diag[i] + shift]; 1259 n = a->i[i+1] - a->i[i]; 1260 PetscLogFlops(2*n-1); 1261 idx = a->j + a->i[i] + shift; 1262 v = a->a + a->i[i] + shift; 1263 sum = b[i]; 1264 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1265 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1266 } 1267 #endif 1268 } 1269 } 1270 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1271 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1272 PetscFunctionReturn(0); 1273 } 1274 1275 #undef __FUNCT__ 1276 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1277 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1278 { 1279 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1280 1281 PetscFunctionBegin; 1282 info->rows_global = (double)A->m; 1283 info->columns_global = (double)A->n; 1284 info->rows_local = (double)A->m; 1285 info->columns_local = (double)A->n; 1286 info->block_size = 1.0; 1287 info->nz_allocated = (double)a->maxnz; 1288 info->nz_used = (double)a->nz; 1289 info->nz_unneeded = (double)(a->maxnz - a->nz); 1290 info->assemblies = (double)A->num_ass; 1291 info->mallocs = (double)a->reallocs; 1292 info->memory = A->mem; 1293 if (A->factor) { 1294 info->fill_ratio_given = A->info.fill_ratio_given; 1295 info->fill_ratio_needed = A->info.fill_ratio_needed; 1296 info->factor_mallocs = A->info.factor_mallocs; 1297 } else { 1298 info->fill_ratio_given = 0; 1299 info->fill_ratio_needed = 0; 1300 info->factor_mallocs = 0; 1301 } 1302 PetscFunctionReturn(0); 1303 } 1304 1305 EXTERN int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,MatLUInfo*,Mat*); 1306 EXTERN int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1307 EXTERN int MatLUFactor_SeqAIJ(Mat,IS,IS,MatLUInfo*); 1308 EXTERN int MatSolve_SeqAIJ(Mat,Vec,Vec); 1309 EXTERN int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1310 EXTERN int MatSolveTranspose_SeqAIJ(Mat,Vec,Vec); 1311 EXTERN int MatSolveTransposeAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1312 1313 #undef __FUNCT__ 1314 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1315 int MatZeroRows_SeqAIJ(Mat A,IS is,PetscScalar *diag) 1316 { 1317 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1318 int i,ierr,N,*rows,m = A->m - 1,shift = a->indexshift; 1319 1320 PetscFunctionBegin; 1321 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1322 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1323 if (a->keepzeroedrows) { 1324 for (i=0; i<N; i++) { 1325 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1326 ierr = PetscMemzero(&a->a[a->i[rows[i]]+shift],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1327 } 1328 if (diag) { 1329 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1330 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1331 for (i=0; i<N; i++) { 1332 a->a[a->diag[rows[i]]] = *diag; 1333 } 1334 } 1335 } else { 1336 if (diag) { 1337 for (i=0; i<N; i++) { 1338 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1339 if (a->ilen[rows[i]] > 0) { 1340 a->ilen[rows[i]] = 1; 1341 a->a[a->i[rows[i]]+shift] = *diag; 1342 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1343 } else { /* in case row was completely empty */ 1344 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1345 } 1346 } 1347 } else { 1348 for (i=0; i<N; i++) { 1349 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1350 a->ilen[rows[i]] = 0; 1351 } 1352 } 1353 } 1354 ierr = ISRestoreIndices(is,&rows);CHKERRQ(ierr); 1355 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1356 PetscFunctionReturn(0); 1357 } 1358 1359 #undef __FUNCT__ 1360 #define __FUNCT__ "MatGetRow_SeqAIJ" 1361 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 1362 { 1363 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1364 int *itmp,i,shift = a->indexshift,ierr; 1365 1366 PetscFunctionBegin; 1367 if (row < 0 || row >= A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %d out of range",row); 1368 1369 *nz = a->i[row+1] - a->i[row]; 1370 if (v) *v = a->a + a->i[row] + shift; 1371 if (idx) { 1372 itmp = a->j + a->i[row] + shift; 1373 if (*nz && shift) { 1374 ierr = PetscMalloc((*nz)*sizeof(int),idx);CHKERRQ(ierr); 1375 for (i=0; i<(*nz); i++) {(*idx)[i] = itmp[i] + shift;} 1376 } else if (*nz) { 1377 *idx = itmp; 1378 } 1379 else *idx = 0; 1380 } 1381 PetscFunctionReturn(0); 1382 } 1383 1384 #undef __FUNCT__ 1385 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1386 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 1387 { 1388 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1389 int ierr; 1390 1391 PetscFunctionBegin; 1392 if (idx) {if (*idx && a->indexshift) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} 1393 PetscFunctionReturn(0); 1394 } 1395 1396 #undef __FUNCT__ 1397 #define __FUNCT__ "MatNorm_SeqAIJ" 1398 int MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1399 { 1400 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1401 PetscScalar *v = a->a; 1402 PetscReal sum = 0.0; 1403 int i,j,shift = a->indexshift,ierr; 1404 1405 PetscFunctionBegin; 1406 if (type == NORM_FROBENIUS) { 1407 for (i=0; i<a->nz; i++) { 1408 #if defined(PETSC_USE_COMPLEX) 1409 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1410 #else 1411 sum += (*v)*(*v); v++; 1412 #endif 1413 } 1414 *nrm = sqrt(sum); 1415 } else if (type == NORM_1) { 1416 PetscReal *tmp; 1417 int *jj = a->j; 1418 ierr = PetscMalloc((A->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1419 ierr = PetscMemzero(tmp,A->n*sizeof(PetscReal));CHKERRQ(ierr); 1420 *nrm = 0.0; 1421 for (j=0; j<a->nz; j++) { 1422 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1423 } 1424 for (j=0; j<A->n; j++) { 1425 if (tmp[j] > *nrm) *nrm = tmp[j]; 1426 } 1427 ierr = PetscFree(tmp);CHKERRQ(ierr); 1428 } else if (type == NORM_INFINITY) { 1429 *nrm = 0.0; 1430 for (j=0; j<A->m; j++) { 1431 v = a->a + a->i[j] + shift; 1432 sum = 0.0; 1433 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1434 sum += PetscAbsScalar(*v); v++; 1435 } 1436 if (sum > *nrm) *nrm = sum; 1437 } 1438 } else { 1439 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1440 } 1441 PetscFunctionReturn(0); 1442 } 1443 1444 #undef __FUNCT__ 1445 #define __FUNCT__ "MatTranspose_SeqAIJ" 1446 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1447 { 1448 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1449 Mat C; 1450 int i,ierr,*aj = a->j,*ai = a->i,m = A->m,len,*col; 1451 int shift = a->indexshift; 1452 PetscScalar *array = a->a; 1453 1454 PetscFunctionBegin; 1455 if (!B && m != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1456 ierr = PetscMalloc((1+A->n)*sizeof(int),&col);CHKERRQ(ierr); 1457 ierr = PetscMemzero(col,(1+A->n)*sizeof(int));CHKERRQ(ierr); 1458 if (shift) { 1459 for (i=0; i<ai[m]-1; i++) aj[i] -= 1; 1460 } 1461 for (i=0; i<ai[m]+shift; i++) col[aj[i]] += 1; 1462 ierr = MatCreateSeqAIJ(A->comm,A->n,m,0,col,&C);CHKERRQ(ierr); 1463 ierr = PetscFree(col);CHKERRQ(ierr); 1464 for (i=0; i<m; i++) { 1465 len = ai[i+1]-ai[i]; 1466 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1467 array += len; 1468 aj += len; 1469 } 1470 if (shift) { 1471 for (i=0; i<ai[m]-1; i++) aj[i] += 1; 1472 } 1473 1474 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1475 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1476 1477 if (B) { 1478 *B = C; 1479 } else { 1480 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1481 } 1482 PetscFunctionReturn(0); 1483 } 1484 1485 #undef __FUNCT__ 1486 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1487 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1488 { 1489 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1490 PetscScalar *l,*r,x,*v; 1491 int ierr,i,j,m = A->m,n = A->n,M,nz = a->nz,*jj,shift = a->indexshift; 1492 1493 PetscFunctionBegin; 1494 if (ll) { 1495 /* The local size is used so that VecMPI can be passed to this routine 1496 by MatDiagonalScale_MPIAIJ */ 1497 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1498 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1499 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1500 v = a->a; 1501 for (i=0; i<m; i++) { 1502 x = l[i]; 1503 M = a->i[i+1] - a->i[i]; 1504 for (j=0; j<M; j++) { (*v++) *= x;} 1505 } 1506 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1507 PetscLogFlops(nz); 1508 } 1509 if (rr) { 1510 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1511 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1512 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1513 v = a->a; jj = a->j; 1514 for (i=0; i<nz; i++) { 1515 (*v++) *= r[*jj++ + shift]; 1516 } 1517 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1518 PetscLogFlops(nz); 1519 } 1520 PetscFunctionReturn(0); 1521 } 1522 1523 #undef __FUNCT__ 1524 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1525 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatReuse scall,Mat *B) 1526 { 1527 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1528 int *smap,i,k,kstart,kend,ierr,oldcols = A->n,*lens; 1529 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1530 int *irow,*icol,nrows,ncols,shift = a->indexshift,*ssmap; 1531 int *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1532 PetscScalar *a_new,*mat_a; 1533 Mat C; 1534 PetscTruth stride; 1535 1536 PetscFunctionBegin; 1537 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1538 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1539 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1540 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1541 1542 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1543 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1544 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1545 1546 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1547 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1548 if (stride && step == 1) { 1549 /* special case of contiguous rows */ 1550 ierr = PetscMalloc((2*nrows+1)*sizeof(int),&lens);CHKERRQ(ierr); 1551 starts = lens + nrows; 1552 /* loop over new rows determining lens and starting points */ 1553 for (i=0; i<nrows; i++) { 1554 kstart = ai[irow[i]]+shift; 1555 kend = kstart + ailen[irow[i]]; 1556 for (k=kstart; k<kend; k++) { 1557 if (aj[k]+shift >= first) { 1558 starts[i] = k; 1559 break; 1560 } 1561 } 1562 sum = 0; 1563 while (k < kend) { 1564 if (aj[k++]+shift >= first+ncols) break; 1565 sum++; 1566 } 1567 lens[i] = sum; 1568 } 1569 /* create submatrix */ 1570 if (scall == MAT_REUSE_MATRIX) { 1571 int n_cols,n_rows; 1572 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1573 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1574 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1575 C = *B; 1576 } else { 1577 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1578 } 1579 c = (Mat_SeqAIJ*)C->data; 1580 1581 /* loop over rows inserting into submatrix */ 1582 a_new = c->a; 1583 j_new = c->j; 1584 i_new = c->i; 1585 i_new[0] = -shift; 1586 for (i=0; i<nrows; i++) { 1587 ii = starts[i]; 1588 lensi = lens[i]; 1589 for (k=0; k<lensi; k++) { 1590 *j_new++ = aj[ii+k] - first; 1591 } 1592 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1593 a_new += lensi; 1594 i_new[i+1] = i_new[i] + lensi; 1595 c->ilen[i] = lensi; 1596 } 1597 ierr = PetscFree(lens);CHKERRQ(ierr); 1598 } else { 1599 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1600 ierr = PetscMalloc((1+oldcols)*sizeof(int),&smap);CHKERRQ(ierr); 1601 ssmap = smap + shift; 1602 ierr = PetscMalloc((1+nrows)*sizeof(int),&lens);CHKERRQ(ierr); 1603 ierr = PetscMemzero(smap,oldcols*sizeof(int));CHKERRQ(ierr); 1604 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1605 /* determine lens of each row */ 1606 for (i=0; i<nrows; i++) { 1607 kstart = ai[irow[i]]+shift; 1608 kend = kstart + a->ilen[irow[i]]; 1609 lens[i] = 0; 1610 for (k=kstart; k<kend; k++) { 1611 if (ssmap[aj[k]]) { 1612 lens[i]++; 1613 } 1614 } 1615 } 1616 /* Create and fill new matrix */ 1617 if (scall == MAT_REUSE_MATRIX) { 1618 PetscTruth equal; 1619 1620 c = (Mat_SeqAIJ *)((*B)->data); 1621 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1622 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(int),&equal);CHKERRQ(ierr); 1623 if (!equal) { 1624 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1625 } 1626 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(int));CHKERRQ(ierr); 1627 C = *B; 1628 } else { 1629 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1630 } 1631 c = (Mat_SeqAIJ *)(C->data); 1632 for (i=0; i<nrows; i++) { 1633 row = irow[i]; 1634 kstart = ai[row]+shift; 1635 kend = kstart + a->ilen[row]; 1636 mat_i = c->i[i]+shift; 1637 mat_j = c->j + mat_i; 1638 mat_a = c->a + mat_i; 1639 mat_ilen = c->ilen + i; 1640 for (k=kstart; k<kend; k++) { 1641 if ((tcol=ssmap[a->j[k]])) { 1642 *mat_j++ = tcol - (!shift); 1643 *mat_a++ = a->a[k]; 1644 (*mat_ilen)++; 1645 1646 } 1647 } 1648 } 1649 /* Free work space */ 1650 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1651 ierr = PetscFree(smap);CHKERRQ(ierr); 1652 ierr = PetscFree(lens);CHKERRQ(ierr); 1653 } 1654 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1655 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1656 1657 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1658 *B = C; 1659 PetscFunctionReturn(0); 1660 } 1661 1662 /* 1663 */ 1664 #undef __FUNCT__ 1665 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1666 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatILUInfo *info) 1667 { 1668 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1669 int ierr; 1670 Mat outA; 1671 PetscTruth row_identity,col_identity; 1672 1673 PetscFunctionBegin; 1674 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1675 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1676 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1677 if (!row_identity || !col_identity) { 1678 SETERRQ(1,"Row and column permutations must be identity for in-place ILU"); 1679 } 1680 1681 outA = inA; 1682 inA->factor = FACTOR_LU; 1683 a->row = row; 1684 a->col = col; 1685 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1686 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1687 1688 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1689 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1690 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1691 PetscLogObjectParent(inA,a->icol); 1692 1693 if (!a->solve_work) { /* this matrix may have been factored before */ 1694 ierr = PetscMalloc((inA->m+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1695 } 1696 1697 if (!a->diag) { 1698 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1699 } 1700 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1701 PetscFunctionReturn(0); 1702 } 1703 1704 #include "petscblaslapack.h" 1705 #undef __FUNCT__ 1706 #define __FUNCT__ "MatScale_SeqAIJ" 1707 int MatScale_SeqAIJ(PetscScalar *alpha,Mat inA) 1708 { 1709 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1710 int one = 1; 1711 1712 PetscFunctionBegin; 1713 BLscal_(&a->nz,alpha,a->a,&one); 1714 PetscLogFlops(a->nz); 1715 PetscFunctionReturn(0); 1716 } 1717 1718 #undef __FUNCT__ 1719 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1720 int MatGetSubMatrices_SeqAIJ(Mat A,int n,IS *irow,IS *icol,MatReuse scall,Mat **B) 1721 { 1722 int ierr,i; 1723 1724 PetscFunctionBegin; 1725 if (scall == MAT_INITIAL_MATRIX) { 1726 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1727 } 1728 1729 for (i=0; i<n; i++) { 1730 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1731 } 1732 PetscFunctionReturn(0); 1733 } 1734 1735 #undef __FUNCT__ 1736 #define __FUNCT__ "MatGetBlockSize_SeqAIJ" 1737 int MatGetBlockSize_SeqAIJ(Mat A,int *bs) 1738 { 1739 PetscFunctionBegin; 1740 *bs = 1; 1741 PetscFunctionReturn(0); 1742 } 1743 1744 #undef __FUNCT__ 1745 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1746 int MatIncreaseOverlap_SeqAIJ(Mat A,int is_max,IS *is,int ov) 1747 { 1748 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1749 int shift,row,i,j,k,l,m,n,*idx,ierr,*nidx,isz,val; 1750 int start,end,*ai,*aj; 1751 PetscBT table; 1752 1753 PetscFunctionBegin; 1754 shift = a->indexshift; 1755 m = A->m; 1756 ai = a->i; 1757 aj = a->j+shift; 1758 1759 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal overlap value used"); 1760 1761 ierr = PetscMalloc((m+1)*sizeof(int),&nidx);CHKERRQ(ierr); 1762 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1763 1764 for (i=0; i<is_max; i++) { 1765 /* Initialize the two local arrays */ 1766 isz = 0; 1767 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1768 1769 /* Extract the indices, assume there can be duplicate entries */ 1770 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1771 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1772 1773 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1774 for (j=0; j<n ; ++j){ 1775 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1776 } 1777 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1778 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1779 1780 k = 0; 1781 for (j=0; j<ov; j++){ /* for each overlap */ 1782 n = isz; 1783 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1784 row = nidx[k]; 1785 start = ai[row]; 1786 end = ai[row+1]; 1787 for (l = start; l<end ; l++){ 1788 val = aj[l] + shift; 1789 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1790 } 1791 } 1792 } 1793 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1794 } 1795 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1796 ierr = PetscFree(nidx);CHKERRQ(ierr); 1797 PetscFunctionReturn(0); 1798 } 1799 1800 /* -------------------------------------------------------------- */ 1801 #undef __FUNCT__ 1802 #define __FUNCT__ "MatPermute_SeqAIJ" 1803 int MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1804 { 1805 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1806 PetscScalar *vwork; 1807 int i,ierr,nz,m = A->m,n = A->n,*cwork; 1808 int *row,*col,*cnew,j,*lens; 1809 IS icolp,irowp; 1810 1811 PetscFunctionBegin; 1812 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1813 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1814 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1815 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1816 1817 /* determine lengths of permuted rows */ 1818 ierr = PetscMalloc((m+1)*sizeof(int),&lens);CHKERRQ(ierr); 1819 for (i=0; i<m; i++) { 1820 lens[row[i]] = a->i[i+1] - a->i[i]; 1821 } 1822 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1823 ierr = PetscFree(lens);CHKERRQ(ierr); 1824 1825 ierr = PetscMalloc(n*sizeof(int),&cnew);CHKERRQ(ierr); 1826 for (i=0; i<m; i++) { 1827 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1828 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1829 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1830 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1831 } 1832 ierr = PetscFree(cnew);CHKERRQ(ierr); 1833 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1834 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1835 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1836 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1837 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1838 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1839 PetscFunctionReturn(0); 1840 } 1841 1842 #undef __FUNCT__ 1843 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1844 int MatPrintHelp_SeqAIJ(Mat A) 1845 { 1846 static PetscTruth called = PETSC_FALSE; 1847 MPI_Comm comm = A->comm; 1848 int ierr; 1849 1850 PetscFunctionBegin; 1851 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1852 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1853 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1854 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1855 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1856 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1857 #if defined(PETSC_HAVE_ESSL) 1858 ierr = (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n");CHKERRQ(ierr); 1859 #endif 1860 #if defined(PETSC_HAVE_LUSOL) 1861 ierr = (*PetscHelpPrintf)(comm," -mat_aij_lusol: Use the Stanford LUSOL sparse factorization and solve.\n");CHKERRQ(ierr); 1862 #endif 1863 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 1864 ierr = (*PetscHelpPrintf)(comm," -mat_aij_matlab: Use Matlab engine sparse LU factorization and solve.\n");CHKERRQ(ierr); 1865 #endif 1866 PetscFunctionReturn(0); 1867 } 1868 EXTERN int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg); 1869 EXTERN int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1870 EXTERN int MatILUDTFactor_SeqAIJ(Mat,MatILUInfo*,IS,IS,Mat*); 1871 #undef __FUNCT__ 1872 #define __FUNCT__ "MatCopy_SeqAIJ" 1873 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1874 { 1875 int ierr; 1876 PetscTruth flg; 1877 1878 PetscFunctionBegin; 1879 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flg);CHKERRQ(ierr); 1880 if (str == SAME_NONZERO_PATTERN && flg) { 1881 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1882 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1883 1884 if (a->i[A->m]+a->indexshift != b->i[B->m]+a->indexshift) { 1885 SETERRQ(1,"Number of nonzeros in two matrices are different"); 1886 } 1887 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m]+a->indexshift)*sizeof(PetscScalar));CHKERRQ(ierr); 1888 } else { 1889 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1890 } 1891 PetscFunctionReturn(0); 1892 } 1893 1894 #undef __FUNCT__ 1895 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1896 int MatSetUpPreallocation_SeqAIJ(Mat A) 1897 { 1898 int ierr; 1899 1900 PetscFunctionBegin; 1901 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1902 PetscFunctionReturn(0); 1903 } 1904 1905 #undef __FUNCT__ 1906 #define __FUNCT__ "MatGetArray_SeqAIJ" 1907 int MatGetArray_SeqAIJ(Mat A,PetscScalar **array) 1908 { 1909 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1910 PetscFunctionBegin; 1911 *array = a->a; 1912 PetscFunctionReturn(0); 1913 } 1914 1915 #undef __FUNCT__ 1916 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1917 int MatRestoreArray_SeqAIJ(Mat A,PetscScalar **array) 1918 { 1919 PetscFunctionBegin; 1920 PetscFunctionReturn(0); 1921 } 1922 1923 #undef __FUNCT__ 1924 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1925 int MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1926 { 1927 int (*f)(void *,Vec,Vec,void*) = (int (*)(void *,Vec,Vec,void *))coloring->f; 1928 int k,ierr,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1929 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1930 PetscScalar *vscale_array; 1931 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1932 Vec w1,w2,w3; 1933 void *fctx = coloring->fctx; 1934 PetscTruth flg; 1935 1936 PetscFunctionBegin; 1937 if (!coloring->w1) { 1938 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1939 PetscLogObjectParent(coloring,coloring->w1); 1940 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1941 PetscLogObjectParent(coloring,coloring->w2); 1942 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1943 PetscLogObjectParent(coloring,coloring->w3); 1944 } 1945 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1946 1947 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1948 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1949 if (flg) { 1950 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 1951 } else { 1952 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1953 } 1954 1955 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1956 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1957 1958 /* 1959 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1960 coloring->F for the coarser grids from the finest 1961 */ 1962 if (coloring->F) { 1963 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1964 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1965 if (m1 != m2) { 1966 coloring->F = 0; 1967 } 1968 } 1969 1970 if (coloring->F) { 1971 w1 = coloring->F; 1972 coloring->F = 0; 1973 } else { 1974 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1975 } 1976 1977 /* 1978 Compute all the scale factors and share with other processors 1979 */ 1980 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 1981 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 1982 for (k=0; k<coloring->ncolors; k++) { 1983 /* 1984 Loop over each column associated with color adding the 1985 perturbation to the vector w3. 1986 */ 1987 for (l=0; l<coloring->ncolumns[k]; l++) { 1988 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1989 dx = xx[col]; 1990 if (dx == 0.0) dx = 1.0; 1991 #if !defined(PETSC_USE_COMPLEX) 1992 if (dx < umin && dx >= 0.0) dx = umin; 1993 else if (dx < 0.0 && dx > -umin) dx = -umin; 1994 #else 1995 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1996 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1997 #endif 1998 dx *= epsilon; 1999 vscale_array[col] = 1.0/dx; 2000 } 2001 } 2002 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2003 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2004 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2005 2006 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2007 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2008 2009 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2010 else vscaleforrow = coloring->columnsforrow; 2011 2012 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2013 /* 2014 Loop over each color 2015 */ 2016 for (k=0; k<coloring->ncolors; k++) { 2017 coloring->currentcolor = k; 2018 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2019 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2020 /* 2021 Loop over each column associated with color adding the 2022 perturbation to the vector w3. 2023 */ 2024 for (l=0; l<coloring->ncolumns[k]; l++) { 2025 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2026 dx = xx[col]; 2027 if (dx == 0.0) dx = 1.0; 2028 #if !defined(PETSC_USE_COMPLEX) 2029 if (dx < umin && dx >= 0.0) dx = umin; 2030 else if (dx < 0.0 && dx > -umin) dx = -umin; 2031 #else 2032 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2033 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2034 #endif 2035 dx *= epsilon; 2036 if (!PetscAbsScalar(dx)) SETERRQ(1,"Computed 0 differencing parameter"); 2037 w3_array[col] += dx; 2038 } 2039 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2040 2041 /* 2042 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2043 */ 2044 2045 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2046 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 2047 2048 /* 2049 Loop over rows of vector, putting results into Jacobian matrix 2050 */ 2051 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2052 for (l=0; l<coloring->nrows[k]; l++) { 2053 row = coloring->rows[k][l]; 2054 col = coloring->columnsforrow[k][l]; 2055 y[row] *= vscale_array[vscaleforrow[k][l]]; 2056 srow = row + start; 2057 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2058 } 2059 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2060 } 2061 coloring->currentcolor = k; 2062 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2063 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2064 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2065 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2066 PetscFunctionReturn(0); 2067 } 2068 2069 #include "petscblaslapack.h" 2070 2071 #undef __FUNCT__ 2072 #define __FUNCT__ "MatAXPY_SeqAIJ" 2073 int MatAXPY_SeqAIJ(PetscScalar *a,Mat X,Mat Y,MatStructure str) 2074 { 2075 int ierr,one=1; 2076 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2077 2078 PetscFunctionBegin; 2079 if (str == SAME_NONZERO_PATTERN) { 2080 BLaxpy_(&x->nz,a,x->a,&one,y->a,&one); 2081 } else { 2082 ierr = MatAXPY_Basic(a,X,Y,str);CHKERRQ(ierr); 2083 } 2084 PetscFunctionReturn(0); 2085 } 2086 2087 2088 /* -------------------------------------------------------------------*/ 2089 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2090 MatGetRow_SeqAIJ, 2091 MatRestoreRow_SeqAIJ, 2092 MatMult_SeqAIJ, 2093 MatMultAdd_SeqAIJ, 2094 MatMultTranspose_SeqAIJ, 2095 MatMultTransposeAdd_SeqAIJ, 2096 MatSolve_SeqAIJ, 2097 MatSolveAdd_SeqAIJ, 2098 MatSolveTranspose_SeqAIJ, 2099 MatSolveTransposeAdd_SeqAIJ, 2100 MatLUFactor_SeqAIJ, 2101 0, 2102 MatRelax_SeqAIJ, 2103 MatTranspose_SeqAIJ, 2104 MatGetInfo_SeqAIJ, 2105 MatEqual_SeqAIJ, 2106 MatGetDiagonal_SeqAIJ, 2107 MatDiagonalScale_SeqAIJ, 2108 MatNorm_SeqAIJ, 2109 0, 2110 MatAssemblyEnd_SeqAIJ, 2111 MatCompress_SeqAIJ, 2112 MatSetOption_SeqAIJ, 2113 MatZeroEntries_SeqAIJ, 2114 MatZeroRows_SeqAIJ, 2115 MatLUFactorSymbolic_SeqAIJ, 2116 MatLUFactorNumeric_SeqAIJ, 2117 0, 2118 MatCholeskyFactorNumeric_SeqAIJ, 2119 MatSetUpPreallocation_SeqAIJ, 2120 MatILUFactorSymbolic_SeqAIJ, 2121 MatICCFactorSymbolic_SeqAIJ, 2122 MatGetArray_SeqAIJ, 2123 MatRestoreArray_SeqAIJ, 2124 MatDuplicate_SeqAIJ, 2125 0, 2126 0, 2127 MatILUFactor_SeqAIJ, 2128 0, 2129 MatAXPY_SeqAIJ, 2130 MatGetSubMatrices_SeqAIJ, 2131 MatIncreaseOverlap_SeqAIJ, 2132 MatGetValues_SeqAIJ, 2133 MatCopy_SeqAIJ, 2134 MatPrintHelp_SeqAIJ, 2135 MatScale_SeqAIJ, 2136 0, 2137 0, 2138 MatILUDTFactor_SeqAIJ, 2139 MatGetBlockSize_SeqAIJ, 2140 MatGetRowIJ_SeqAIJ, 2141 MatRestoreRowIJ_SeqAIJ, 2142 MatGetColumnIJ_SeqAIJ, 2143 MatRestoreColumnIJ_SeqAIJ, 2144 MatFDColoringCreate_SeqAIJ, 2145 0, 2146 0, 2147 MatPermute_SeqAIJ, 2148 0, 2149 0, 2150 MatDestroy_SeqAIJ, 2151 MatView_SeqAIJ, 2152 MatGetPetscMaps_Petsc, 2153 0, 2154 0, 2155 0, 2156 0, 2157 0, 2158 0, 2159 0, 2160 0, 2161 MatSetColoring_SeqAIJ, 2162 MatSetValuesAdic_SeqAIJ, 2163 MatSetValuesAdifor_SeqAIJ, 2164 MatFDColoringApply_SeqAIJ}; 2165 2166 EXTERN_C_BEGIN 2167 #undef __FUNCT__ 2168 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2169 2170 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 2171 { 2172 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2173 int i,nz,n; 2174 2175 PetscFunctionBegin; 2176 2177 nz = aij->maxnz; 2178 n = mat->n; 2179 for (i=0; i<nz; i++) { 2180 aij->j[i] = indices[i]; 2181 } 2182 aij->nz = nz; 2183 for (i=0; i<n; i++) { 2184 aij->ilen[i] = aij->imax[i]; 2185 } 2186 2187 PetscFunctionReturn(0); 2188 } 2189 EXTERN_C_END 2190 2191 #undef __FUNCT__ 2192 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2193 /*@ 2194 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2195 in the matrix. 2196 2197 Input Parameters: 2198 + mat - the SeqAIJ matrix 2199 - indices - the column indices 2200 2201 Level: advanced 2202 2203 Notes: 2204 This can be called if you have precomputed the nonzero structure of the 2205 matrix and want to provide it to the matrix object to improve the performance 2206 of the MatSetValues() operation. 2207 2208 You MUST have set the correct numbers of nonzeros per row in the call to 2209 MatCreateSeqAIJ(). 2210 2211 MUST be called before any calls to MatSetValues(); 2212 2213 The indices should start with zero, not one. 2214 2215 @*/ 2216 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 2217 { 2218 int ierr,(*f)(Mat,int *); 2219 2220 PetscFunctionBegin; 2221 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2222 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2223 if (f) { 2224 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2225 } else { 2226 SETERRQ(1,"Wrong type of matrix to set column indices"); 2227 } 2228 PetscFunctionReturn(0); 2229 } 2230 2231 /* ----------------------------------------------------------------------------------------*/ 2232 2233 EXTERN_C_BEGIN 2234 #undef __FUNCT__ 2235 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2236 int MatStoreValues_SeqAIJ(Mat mat) 2237 { 2238 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2239 size_t nz = aij->i[mat->m]+aij->indexshift,ierr; 2240 2241 PetscFunctionBegin; 2242 if (aij->nonew != 1) { 2243 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2244 } 2245 2246 /* allocate space for values if not already there */ 2247 if (!aij->saved_values) { 2248 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2249 } 2250 2251 /* copy values over */ 2252 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2253 PetscFunctionReturn(0); 2254 } 2255 EXTERN_C_END 2256 2257 #undef __FUNCT__ 2258 #define __FUNCT__ "MatStoreValues" 2259 /*@ 2260 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2261 example, reuse of the linear part of a Jacobian, while recomputing the 2262 nonlinear portion. 2263 2264 Collect on Mat 2265 2266 Input Parameters: 2267 . mat - the matrix (currently on AIJ matrices support this option) 2268 2269 Level: advanced 2270 2271 Common Usage, with SNESSolve(): 2272 $ Create Jacobian matrix 2273 $ Set linear terms into matrix 2274 $ Apply boundary conditions to matrix, at this time matrix must have 2275 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2276 $ boundary conditions again will not change the nonzero structure 2277 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2278 $ ierr = MatStoreValues(mat); 2279 $ Call SNESSetJacobian() with matrix 2280 $ In your Jacobian routine 2281 $ ierr = MatRetrieveValues(mat); 2282 $ Set nonlinear terms in matrix 2283 2284 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2285 $ // build linear portion of Jacobian 2286 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2287 $ ierr = MatStoreValues(mat); 2288 $ loop over nonlinear iterations 2289 $ ierr = MatRetrieveValues(mat); 2290 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2291 $ // call MatAssemblyBegin/End() on matrix 2292 $ Solve linear system with Jacobian 2293 $ endloop 2294 2295 Notes: 2296 Matrix must already be assemblied before calling this routine 2297 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2298 calling this routine. 2299 2300 .seealso: MatRetrieveValues() 2301 2302 @*/ 2303 int MatStoreValues(Mat mat) 2304 { 2305 int ierr,(*f)(Mat); 2306 2307 PetscFunctionBegin; 2308 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2309 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2310 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2311 2312 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2313 if (f) { 2314 ierr = (*f)(mat);CHKERRQ(ierr); 2315 } else { 2316 SETERRQ(1,"Wrong type of matrix to store values"); 2317 } 2318 PetscFunctionReturn(0); 2319 } 2320 2321 EXTERN_C_BEGIN 2322 #undef __FUNCT__ 2323 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2324 int MatRetrieveValues_SeqAIJ(Mat mat) 2325 { 2326 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2327 int nz = aij->i[mat->m]+aij->indexshift,ierr; 2328 2329 PetscFunctionBegin; 2330 if (aij->nonew != 1) { 2331 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2332 } 2333 if (!aij->saved_values) { 2334 SETERRQ(1,"Must call MatStoreValues(A);first"); 2335 } 2336 2337 /* copy values over */ 2338 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2339 PetscFunctionReturn(0); 2340 } 2341 EXTERN_C_END 2342 2343 #undef __FUNCT__ 2344 #define __FUNCT__ "MatRetrieveValues" 2345 /*@ 2346 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2347 example, reuse of the linear part of a Jacobian, while recomputing the 2348 nonlinear portion. 2349 2350 Collect on Mat 2351 2352 Input Parameters: 2353 . mat - the matrix (currently on AIJ matrices support this option) 2354 2355 Level: advanced 2356 2357 .seealso: MatStoreValues() 2358 2359 @*/ 2360 int MatRetrieveValues(Mat mat) 2361 { 2362 int ierr,(*f)(Mat); 2363 2364 PetscFunctionBegin; 2365 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2366 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2367 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2368 2369 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2370 if (f) { 2371 ierr = (*f)(mat);CHKERRQ(ierr); 2372 } else { 2373 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2374 } 2375 PetscFunctionReturn(0); 2376 } 2377 2378 /* 2379 This allows SeqAIJ matrices to be passed to the matlab engine 2380 */ 2381 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2382 #include "engine.h" /* Matlab include file */ 2383 #include "mex.h" /* Matlab include file */ 2384 EXTERN_C_BEGIN 2385 #undef __FUNCT__ 2386 #define __FUNCT__ "MatMatlabEnginePut_SeqAIJ" 2387 int MatMatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine) 2388 { 2389 int ierr,i,*ai,*aj; 2390 Mat B = (Mat)obj; 2391 mxArray *mat; 2392 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2393 2394 PetscFunctionBegin; 2395 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2396 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2397 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2398 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2399 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2400 2401 /* Matlab indices start at 0 for sparse (what a surprise) */ 2402 if (aij->indexshift) { 2403 ai = mxGetJc(mat); 2404 for (i=0; i<B->m+1; i++) { 2405 ai[i]--; 2406 } 2407 aj = mxGetIr(mat); 2408 for (i=0; i<aij->nz; i++) { 2409 aj[i]--; 2410 } 2411 } 2412 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2413 mxSetName(mat,obj->name); 2414 engPutArray((Engine *)mengine,mat); 2415 PetscFunctionReturn(0); 2416 } 2417 EXTERN_C_END 2418 2419 EXTERN_C_BEGIN 2420 #undef __FUNCT__ 2421 #define __FUNCT__ "MatMatlabEngineGet_SeqAIJ" 2422 int MatMatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine) 2423 { 2424 int ierr,ii; 2425 Mat mat = (Mat)obj; 2426 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2427 mxArray *mmat; 2428 2429 PetscFunctionBegin; 2430 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2431 aij->indexshift = 0; 2432 2433 mmat = engGetArray((Engine *)mengine,obj->name); 2434 2435 aij->nz = (mxGetJc(mmat))[mat->m]; 2436 ierr = PetscMalloc(((size_t) aij->nz)*(sizeof(int)+sizeof(PetscScalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 2437 aij->j = (int*)(aij->a + aij->nz); 2438 aij->i = aij->j + aij->nz; 2439 aij->singlemalloc = PETSC_TRUE; 2440 aij->freedata = PETSC_TRUE; 2441 2442 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2443 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2444 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 2445 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 2446 2447 for (ii=0; ii<mat->m; ii++) { 2448 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 2449 } 2450 2451 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2452 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2453 2454 PetscFunctionReturn(0); 2455 } 2456 EXTERN_C_END 2457 #endif 2458 2459 /* --------------------------------------------------------------------------------*/ 2460 #undef __FUNCT__ 2461 #define __FUNCT__ "MatCreateSeqAIJ" 2462 /*@C 2463 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2464 (the default parallel PETSc format). For good matrix assembly performance 2465 the user should preallocate the matrix storage by setting the parameter nz 2466 (or the array nnz). By setting these parameters accurately, performance 2467 during matrix assembly can be increased by more than a factor of 50. 2468 2469 Collective on MPI_Comm 2470 2471 Input Parameters: 2472 + comm - MPI communicator, set to PETSC_COMM_SELF 2473 . m - number of rows 2474 . n - number of columns 2475 . nz - number of nonzeros per row (same for all rows) 2476 - nnz - array containing the number of nonzeros in the various rows 2477 (possibly different for each row) or PETSC_NULL 2478 2479 Output Parameter: 2480 . A - the matrix 2481 2482 Notes: 2483 The AIJ format (also called the Yale sparse matrix format or 2484 compressed row storage), is fully compatible with standard Fortran 77 2485 storage. That is, the stored row and column indices can begin at 2486 either one (as in Fortran) or zero. See the users' manual for details. 2487 2488 Specify the preallocated storage with either nz or nnz (not both). 2489 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2490 allocation. For large problems you MUST preallocate memory or you 2491 will get TERRIBLE performance, see the users' manual chapter on matrices. 2492 2493 By default, this format uses inodes (identical nodes) when possible, to 2494 improve numerical efficiency of matrix-vector products and solves. We 2495 search for consecutive rows with the same nonzero structure, thereby 2496 reusing matrix information to achieve increased efficiency. 2497 2498 Options Database Keys: 2499 + -mat_aij_no_inode - Do not use inodes 2500 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2501 - -mat_aij_oneindex - Internally use indexing starting at 1 2502 rather than 0. Note that when calling MatSetValues(), 2503 the user still MUST index entries starting at 0! 2504 2505 Level: intermediate 2506 2507 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2508 2509 @*/ 2510 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz,Mat *A) 2511 { 2512 int ierr; 2513 2514 PetscFunctionBegin; 2515 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2516 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2517 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2518 PetscFunctionReturn(0); 2519 } 2520 2521 #define SKIP_ALLOCATION -4 2522 2523 #undef __FUNCT__ 2524 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2525 /*@C 2526 MatSeqAIJSetPreallocation - For good matrix assembly performance 2527 the user should preallocate the matrix storage by setting the parameter nz 2528 (or the array nnz). By setting these parameters accurately, performance 2529 during matrix assembly can be increased by more than a factor of 50. 2530 2531 Collective on MPI_Comm 2532 2533 Input Parameters: 2534 + comm - MPI communicator, set to PETSC_COMM_SELF 2535 . m - number of rows 2536 . n - number of columns 2537 . nz - number of nonzeros per row (same for all rows) 2538 - nnz - array containing the number of nonzeros in the various rows 2539 (possibly different for each row) or PETSC_NULL 2540 2541 Output Parameter: 2542 . A - the matrix 2543 2544 Notes: 2545 The AIJ format (also called the Yale sparse matrix format or 2546 compressed row storage), is fully compatible with standard Fortran 77 2547 storage. That is, the stored row and column indices can begin at 2548 either one (as in Fortran) or zero. See the users' manual for details. 2549 2550 Specify the preallocated storage with either nz or nnz (not both). 2551 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2552 allocation. For large problems you MUST preallocate memory or you 2553 will get TERRIBLE performance, see the users' manual chapter on matrices. 2554 2555 By default, this format uses inodes (identical nodes) when possible, to 2556 improve numerical efficiency of matrix-vector products and solves. We 2557 search for consecutive rows with the same nonzero structure, thereby 2558 reusing matrix information to achieve increased efficiency. 2559 2560 Options Database Keys: 2561 + -mat_aij_no_inode - Do not use inodes 2562 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2563 - -mat_aij_oneindex - Internally use indexing starting at 1 2564 rather than 0. Note that when calling MatSetValues(), 2565 the user still MUST index entries starting at 0! 2566 2567 Level: intermediate 2568 2569 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2570 2571 @*/ 2572 int MatSeqAIJSetPreallocation(Mat B,int nz,int *nnz) 2573 { 2574 Mat_SeqAIJ *b; 2575 size_t len = 0; 2576 PetscTruth flg2,skipallocation = PETSC_FALSE; 2577 int i,ierr; 2578 2579 PetscFunctionBegin; 2580 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flg2);CHKERRQ(ierr); 2581 if (!flg2) PetscFunctionReturn(0); 2582 2583 if (nz == SKIP_ALLOCATION) { 2584 skipallocation = PETSC_TRUE; 2585 nz = 0; 2586 } 2587 2588 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2589 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2590 if (nnz) { 2591 for (i=0; i<B->m; i++) { 2592 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2593 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); 2594 } 2595 } 2596 2597 B->preallocated = PETSC_TRUE; 2598 b = (Mat_SeqAIJ*)B->data; 2599 2600 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2601 if (!nnz) { 2602 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2603 else if (nz <= 0) nz = 1; 2604 for (i=0; i<B->m; i++) b->imax[i] = nz; 2605 nz = nz*B->m; 2606 } else { 2607 nz = 0; 2608 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2609 } 2610 2611 if (!skipallocation) { 2612 /* allocate the matrix space */ 2613 len = ((size_t) nz)*(sizeof(int) + sizeof(PetscScalar)) + (B->m+1)*sizeof(int); 2614 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2615 b->j = (int*)(b->a + nz); 2616 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2617 b->i = b->j + nz; 2618 b->i[0] = -b->indexshift; 2619 for (i=1; i<B->m+1; i++) { 2620 b->i[i] = b->i[i-1] + b->imax[i-1]; 2621 } 2622 b->singlemalloc = PETSC_TRUE; 2623 b->freedata = PETSC_TRUE; 2624 } else { 2625 b->freedata = PETSC_FALSE; 2626 } 2627 2628 /* b->ilen will count nonzeros in each row so far. */ 2629 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2630 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2631 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2632 2633 b->nz = 0; 2634 b->maxnz = nz; 2635 B->info.nz_unneeded = (double)b->maxnz; 2636 PetscFunctionReturn(0); 2637 } 2638 2639 EXTERN int RegisterApplyPtAPRoutines_Private(Mat); 2640 2641 EXTERN_C_BEGIN 2642 extern int MatConvert_SeqAIJ_SeqSBAIJ(Mat,MatType,Mat*); 2643 extern int MatConvert_SeqAIJ_SeqBAIJ(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 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2730 "MatConvert_SeqAIJ_SeqBAIJ", 2731 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2732 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2733 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2734 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2735 #endif 2736 ierr = RegisterApplyPtAPRoutines_Private(B);CHKERRQ(ierr); 2737 PetscFunctionReturn(0); 2738 } 2739 EXTERN_C_END 2740 2741 #undef __FUNCT__ 2742 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2743 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2744 { 2745 Mat C; 2746 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2747 int i,m = A->m,shift = a->indexshift,ierr; 2748 size_t len; 2749 2750 PetscFunctionBegin; 2751 *B = 0; 2752 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2753 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 2754 c = (Mat_SeqAIJ*)C->data; 2755 2756 C->factor = A->factor; 2757 c->row = 0; 2758 c->col = 0; 2759 c->icol = 0; 2760 c->indexshift = shift; 2761 c->keepzeroedrows = a->keepzeroedrows; 2762 C->assembled = PETSC_TRUE; 2763 2764 C->M = A->m; 2765 C->N = A->n; 2766 2767 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2768 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2769 for (i=0; i<m; i++) { 2770 c->imax[i] = a->imax[i]; 2771 c->ilen[i] = a->ilen[i]; 2772 } 2773 2774 /* allocate the matrix space */ 2775 c->singlemalloc = PETSC_TRUE; 2776 len = ((size_t) (m+1))*sizeof(int)+(a->i[m])*(sizeof(PetscScalar)+sizeof(int)); 2777 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2778 c->j = (int*)(c->a + a->i[m] + shift); 2779 c->i = c->j + a->i[m] + shift; 2780 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2781 if (m > 0) { 2782 ierr = PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int));CHKERRQ(ierr); 2783 if (cpvalues == MAT_COPY_VALUES) { 2784 ierr = PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(PetscScalar));CHKERRQ(ierr); 2785 } else { 2786 ierr = PetscMemzero(c->a,(a->i[m]+shift)*sizeof(PetscScalar));CHKERRQ(ierr); 2787 } 2788 } 2789 2790 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2791 c->sorted = a->sorted; 2792 c->roworiented = a->roworiented; 2793 c->nonew = a->nonew; 2794 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2795 c->saved_values = 0; 2796 c->idiag = 0; 2797 c->ssor = 0; 2798 c->ignorezeroentries = a->ignorezeroentries; 2799 c->freedata = PETSC_TRUE; 2800 2801 if (a->diag) { 2802 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2803 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2804 for (i=0; i<m; i++) { 2805 c->diag[i] = a->diag[i]; 2806 } 2807 } else c->diag = 0; 2808 c->inode.use = a->inode.use; 2809 c->inode.limit = a->inode.limit; 2810 c->inode.max_limit = a->inode.max_limit; 2811 if (a->inode.size){ 2812 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2813 c->inode.node_count = a->inode.node_count; 2814 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2815 } else { 2816 c->inode.size = 0; 2817 c->inode.node_count = 0; 2818 } 2819 c->nz = a->nz; 2820 c->maxnz = a->maxnz; 2821 c->solve_work = 0; 2822 C->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2823 C->preallocated = PETSC_TRUE; 2824 2825 *B = C; 2826 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2827 PetscFunctionReturn(0); 2828 } 2829 2830 EXTERN_C_BEGIN 2831 #undef __FUNCT__ 2832 #define __FUNCT__ "MatLoad_SeqAIJ" 2833 int MatLoad_SeqAIJ(PetscViewer viewer,MatType type,Mat *A) 2834 { 2835 Mat_SeqAIJ *a; 2836 Mat B; 2837 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N,shift; 2838 MPI_Comm comm; 2839 2840 PetscFunctionBegin; 2841 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2842 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2843 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2844 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2845 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2846 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2847 M = header[1]; N = header[2]; nz = header[3]; 2848 2849 if (nz < 0) { 2850 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2851 } 2852 2853 /* read in row lengths */ 2854 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2855 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2856 2857 /* create our matrix */ 2858 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A);CHKERRQ(ierr); 2859 B = *A; 2860 a = (Mat_SeqAIJ*)B->data; 2861 shift = a->indexshift; 2862 2863 /* read in column indices and adjust for Fortran indexing*/ 2864 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2865 if (shift) { 2866 for (i=0; i<nz; i++) { 2867 a->j[i] += 1; 2868 } 2869 } 2870 2871 /* read in nonzero values */ 2872 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2873 2874 /* set matrix "i" values */ 2875 a->i[0] = -shift; 2876 for (i=1; i<= M; i++) { 2877 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2878 a->ilen[i-1] = rowlengths[i-1]; 2879 } 2880 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2881 2882 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2883 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2884 PetscFunctionReturn(0); 2885 } 2886 EXTERN_C_END 2887 2888 #undef __FUNCT__ 2889 #define __FUNCT__ "MatEqual_SeqAIJ" 2890 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2891 { 2892 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2893 int ierr; 2894 PetscTruth flag; 2895 2896 PetscFunctionBegin; 2897 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flag);CHKERRQ(ierr); 2898 if (!flag) SETERRQ(PETSC_ERR_ARG_INCOMP,"Matrices must be same type"); 2899 2900 /* If the matrix dimensions are not equal,or no of nonzeros or shift */ 2901 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)|| (a->indexshift != b->indexshift)) { 2902 *flg = PETSC_FALSE; 2903 PetscFunctionReturn(0); 2904 } 2905 2906 /* if the a->i are the same */ 2907 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2908 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2909 2910 /* if a->j are the same */ 2911 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2912 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2913 2914 /* if a->a are the same */ 2915 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2916 2917 PetscFunctionReturn(0); 2918 2919 } 2920 2921 #undef __FUNCT__ 2922 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2923 /*@C 2924 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2925 provided by the user. 2926 2927 Coolective on MPI_Comm 2928 2929 Input Parameters: 2930 + comm - must be an MPI communicator of size 1 2931 . m - number of rows 2932 . n - number of columns 2933 . i - row indices 2934 . j - column indices 2935 - a - matrix values 2936 2937 Output Parameter: 2938 . mat - the matrix 2939 2940 Level: intermediate 2941 2942 Notes: 2943 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2944 once the matrix is destroyed 2945 2946 You cannot set new nonzero locations into this matrix, that will generate an error. 2947 2948 The i and j indices can be either 0- or 1 based 2949 2950 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2951 2952 @*/ 2953 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,PetscScalar *a,Mat *mat) 2954 { 2955 int ierr,ii; 2956 Mat_SeqAIJ *aij; 2957 2958 PetscFunctionBegin; 2959 ierr = MatCreateSeqAIJ(comm,m,n,SKIP_ALLOCATION,0,mat);CHKERRQ(ierr); 2960 aij = (Mat_SeqAIJ*)(*mat)->data; 2961 2962 if (i[0] == 1) { 2963 aij->indexshift = -1; 2964 } else if (i[0]) { 2965 SETERRQ(1,"i (row indices) do not start with 0 or 1"); 2966 } 2967 aij->i = i; 2968 aij->j = j; 2969 aij->a = a; 2970 aij->singlemalloc = PETSC_FALSE; 2971 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2972 aij->freedata = PETSC_FALSE; 2973 2974 for (ii=0; ii<m; ii++) { 2975 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2976 #if defined(PETSC_USE_BOPT_g) 2977 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]); 2978 #endif 2979 } 2980 #if defined(PETSC_USE_BOPT_g) 2981 for (ii=0; ii<aij->i[m]; ii++) { 2982 if (j[ii] < -aij->indexshift) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 2983 if (j[ii] > n - 1 -aij->indexshift) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 2984 } 2985 #endif 2986 2987 /* changes indices to start at 0 */ 2988 if (i[0]) { 2989 aij->indexshift = 0; 2990 for (ii=0; ii<m; ii++) { 2991 i[ii]--; 2992 } 2993 for (ii=0; ii<i[m]; ii++) { 2994 j[ii]--; 2995 } 2996 } 2997 2998 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2999 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3000 PetscFunctionReturn(0); 3001 } 3002 3003 #undef __FUNCT__ 3004 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3005 int MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3006 { 3007 int ierr; 3008 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3009 3010 PetscFunctionBegin; 3011 if (coloring->ctype == IS_COLORING_LOCAL) { 3012 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3013 a->coloring = coloring; 3014 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3015 int *colors,i,*larray; 3016 ISColoring ocoloring; 3017 3018 /* set coloring for diagonal portion */ 3019 ierr = PetscMalloc((A->n+1)*sizeof(int),&larray);CHKERRQ(ierr); 3020 for (i=0; i<A->n; i++) { 3021 larray[i] = i; 3022 } 3023 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3024 ierr = PetscMalloc((A->n+1)*sizeof(int),&colors);CHKERRQ(ierr); 3025 for (i=0; i<A->n; i++) { 3026 colors[i] = coloring->colors[larray[i]]; 3027 } 3028 ierr = PetscFree(larray);CHKERRQ(ierr); 3029 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 3030 a->coloring = ocoloring; 3031 } 3032 PetscFunctionReturn(0); 3033 } 3034 3035 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 3036 EXTERN_C_BEGIN 3037 #include "adic/ad_utils.h" 3038 EXTERN_C_END 3039 3040 #undef __FUNCT__ 3041 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3042 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3043 { 3044 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3045 int m = A->m,*ii = a->i,*jj = a->j,nz,i,*color,j,nlen; 3046 PetscScalar *v = a->a,*values; 3047 char *cadvalues = (char *)advalues; 3048 3049 PetscFunctionBegin; 3050 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3051 nlen = PetscADGetDerivTypeSize(); 3052 color = a->coloring->colors; 3053 /* loop over rows */ 3054 for (i=0; i<m; i++) { 3055 nz = ii[i+1] - ii[i]; 3056 /* loop over columns putting computed value into matrix */ 3057 values = PetscADGetGradArray(cadvalues); 3058 for (j=0; j<nz; j++) { 3059 *v++ = values[color[*jj++]]; 3060 } 3061 cadvalues += nlen; /* jump to next row of derivatives */ 3062 } 3063 PetscFunctionReturn(0); 3064 } 3065 3066 #else 3067 3068 #undef __FUNCT__ 3069 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3070 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3071 { 3072 PetscFunctionBegin; 3073 SETERRQ(1,"PETSc installed without ADIC"); 3074 } 3075 3076 #endif 3077 3078 #undef __FUNCT__ 3079 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3080 int MatSetValuesAdifor_SeqAIJ(Mat A,int nl,void *advalues) 3081 { 3082 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3083 int m = A->m,*ii = a->i,*jj = a->j,nz,i,*color,j; 3084 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3085 3086 PetscFunctionBegin; 3087 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3088 color = a->coloring->colors; 3089 /* loop over rows */ 3090 for (i=0; i<m; i++) { 3091 nz = ii[i+1] - ii[i]; 3092 /* loop over columns putting computed value into matrix */ 3093 for (j=0; j<nz; j++) { 3094 *v++ = values[color[*jj++]]; 3095 } 3096 values += nl; /* jump to next row of derivatives */ 3097 } 3098 PetscFunctionReturn(0); 3099 } 3100 3101