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