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