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