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