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 2050 #undef __FUNCT__ 2051 #define __FUNCT__ "MatAXPY_SeqAIJ" 2052 int MatAXPY_SeqAIJ(PetscScalar *a,Mat X,Mat Y,MatStructure str) 2053 { 2054 int ierr,one=1; 2055 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2056 2057 PetscFunctionBegin; 2058 if (str == SAME_NONZERO_PATTERN) { 2059 BLaxpy_(&x->nz,a,x->a,&one,y->a,&one); 2060 } else { 2061 ierr = MatAXPY_Basic(a,X,Y,str);CHKERRQ(ierr); 2062 } 2063 PetscFunctionReturn(0); 2064 } 2065 2066 2067 /* -------------------------------------------------------------------*/ 2068 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2069 MatGetRow_SeqAIJ, 2070 MatRestoreRow_SeqAIJ, 2071 MatMult_SeqAIJ, 2072 MatMultAdd_SeqAIJ, 2073 MatMultTranspose_SeqAIJ, 2074 MatMultTransposeAdd_SeqAIJ, 2075 MatSolve_SeqAIJ, 2076 MatSolveAdd_SeqAIJ, 2077 MatSolveTranspose_SeqAIJ, 2078 MatSolveTransposeAdd_SeqAIJ, 2079 MatLUFactor_SeqAIJ, 2080 0, 2081 MatRelax_SeqAIJ, 2082 MatTranspose_SeqAIJ, 2083 MatGetInfo_SeqAIJ, 2084 MatEqual_SeqAIJ, 2085 MatGetDiagonal_SeqAIJ, 2086 MatDiagonalScale_SeqAIJ, 2087 MatNorm_SeqAIJ, 2088 0, 2089 MatAssemblyEnd_SeqAIJ, 2090 MatCompress_SeqAIJ, 2091 MatSetOption_SeqAIJ, 2092 MatZeroEntries_SeqAIJ, 2093 MatZeroRows_SeqAIJ, 2094 MatLUFactorSymbolic_SeqAIJ, 2095 MatLUFactorNumeric_SeqAIJ, 2096 0, 2097 MatCholeskyFactorNumeric_SeqAIJ, 2098 MatSetUpPreallocation_SeqAIJ, 2099 MatILUFactorSymbolic_SeqAIJ, 2100 MatICCFactorSymbolic_SeqAIJ, 2101 MatGetArray_SeqAIJ, 2102 MatRestoreArray_SeqAIJ, 2103 MatDuplicate_SeqAIJ, 2104 0, 2105 0, 2106 MatILUFactor_SeqAIJ, 2107 0, 2108 MatAXPY_SeqAIJ, 2109 MatGetSubMatrices_SeqAIJ, 2110 MatIncreaseOverlap_SeqAIJ, 2111 MatGetValues_SeqAIJ, 2112 MatCopy_SeqAIJ, 2113 MatPrintHelp_SeqAIJ, 2114 MatScale_SeqAIJ, 2115 0, 2116 0, 2117 MatILUDTFactor_SeqAIJ, 2118 MatGetBlockSize_SeqAIJ, 2119 MatGetRowIJ_SeqAIJ, 2120 MatRestoreRowIJ_SeqAIJ, 2121 MatGetColumnIJ_SeqAIJ, 2122 MatRestoreColumnIJ_SeqAIJ, 2123 MatFDColoringCreate_SeqAIJ, 2124 0, 2125 0, 2126 MatPermute_SeqAIJ, 2127 0, 2128 0, 2129 MatDestroy_SeqAIJ, 2130 MatView_SeqAIJ, 2131 MatGetPetscMaps_Petsc, 2132 0, 2133 0, 2134 0, 2135 0, 2136 0, 2137 0, 2138 0, 2139 0, 2140 MatSetColoring_SeqAIJ, 2141 MatSetValuesAdic_SeqAIJ, 2142 MatSetValuesAdifor_SeqAIJ, 2143 MatFDColoringApply_SeqAIJ}; 2144 2145 EXTERN_C_BEGIN 2146 #undef __FUNCT__ 2147 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2148 2149 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 2150 { 2151 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2152 int i,nz,n; 2153 2154 PetscFunctionBegin; 2155 2156 nz = aij->maxnz; 2157 n = mat->n; 2158 for (i=0; i<nz; i++) { 2159 aij->j[i] = indices[i]; 2160 } 2161 aij->nz = nz; 2162 for (i=0; i<n; i++) { 2163 aij->ilen[i] = aij->imax[i]; 2164 } 2165 2166 PetscFunctionReturn(0); 2167 } 2168 EXTERN_C_END 2169 2170 #undef __FUNCT__ 2171 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2172 /*@ 2173 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2174 in the matrix. 2175 2176 Input Parameters: 2177 + mat - the SeqAIJ matrix 2178 - indices - the column indices 2179 2180 Level: advanced 2181 2182 Notes: 2183 This can be called if you have precomputed the nonzero structure of the 2184 matrix and want to provide it to the matrix object to improve the performance 2185 of the MatSetValues() operation. 2186 2187 You MUST have set the correct numbers of nonzeros per row in the call to 2188 MatCreateSeqAIJ(). 2189 2190 MUST be called before any calls to MatSetValues(); 2191 2192 The indices should start with zero, not one. 2193 2194 @*/ 2195 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 2196 { 2197 int ierr,(*f)(Mat,int *); 2198 2199 PetscFunctionBegin; 2200 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2201 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2202 if (f) { 2203 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2204 } else { 2205 SETERRQ(1,"Wrong type of matrix to set column indices"); 2206 } 2207 PetscFunctionReturn(0); 2208 } 2209 2210 /* ----------------------------------------------------------------------------------------*/ 2211 2212 EXTERN_C_BEGIN 2213 #undef __FUNCT__ 2214 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2215 int MatStoreValues_SeqAIJ(Mat mat) 2216 { 2217 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2218 size_t nz = aij->i[mat->m]+aij->indexshift,ierr; 2219 2220 PetscFunctionBegin; 2221 if (aij->nonew != 1) { 2222 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2223 } 2224 2225 /* allocate space for values if not already there */ 2226 if (!aij->saved_values) { 2227 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2228 } 2229 2230 /* copy values over */ 2231 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2232 PetscFunctionReturn(0); 2233 } 2234 EXTERN_C_END 2235 2236 #undef __FUNCT__ 2237 #define __FUNCT__ "MatStoreValues" 2238 /*@ 2239 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2240 example, reuse of the linear part of a Jacobian, while recomputing the 2241 nonlinear portion. 2242 2243 Collect on Mat 2244 2245 Input Parameters: 2246 . mat - the matrix (currently on AIJ matrices support this option) 2247 2248 Level: advanced 2249 2250 Common Usage, with SNESSolve(): 2251 $ Create Jacobian matrix 2252 $ Set linear terms into matrix 2253 $ Apply boundary conditions to matrix, at this time matrix must have 2254 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2255 $ boundary conditions again will not change the nonzero structure 2256 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2257 $ ierr = MatStoreValues(mat); 2258 $ Call SNESSetJacobian() with matrix 2259 $ In your Jacobian routine 2260 $ ierr = MatRetrieveValues(mat); 2261 $ Set nonlinear terms in matrix 2262 2263 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2264 $ // build linear portion of Jacobian 2265 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2266 $ ierr = MatStoreValues(mat); 2267 $ loop over nonlinear iterations 2268 $ ierr = MatRetrieveValues(mat); 2269 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2270 $ // call MatAssemblyBegin/End() on matrix 2271 $ Solve linear system with Jacobian 2272 $ endloop 2273 2274 Notes: 2275 Matrix must already be assemblied before calling this routine 2276 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2277 calling this routine. 2278 2279 .seealso: MatRetrieveValues() 2280 2281 @*/ 2282 int MatStoreValues(Mat mat) 2283 { 2284 int ierr,(*f)(Mat); 2285 2286 PetscFunctionBegin; 2287 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2288 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2289 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2290 2291 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2292 if (f) { 2293 ierr = (*f)(mat);CHKERRQ(ierr); 2294 } else { 2295 SETERRQ(1,"Wrong type of matrix to store values"); 2296 } 2297 PetscFunctionReturn(0); 2298 } 2299 2300 EXTERN_C_BEGIN 2301 #undef __FUNCT__ 2302 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2303 int MatRetrieveValues_SeqAIJ(Mat mat) 2304 { 2305 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2306 int nz = aij->i[mat->m]+aij->indexshift,ierr; 2307 2308 PetscFunctionBegin; 2309 if (aij->nonew != 1) { 2310 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2311 } 2312 if (!aij->saved_values) { 2313 SETERRQ(1,"Must call MatStoreValues(A);first"); 2314 } 2315 2316 /* copy values over */ 2317 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2318 PetscFunctionReturn(0); 2319 } 2320 EXTERN_C_END 2321 2322 #undef __FUNCT__ 2323 #define __FUNCT__ "MatRetrieveValues" 2324 /*@ 2325 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2326 example, reuse of the linear part of a Jacobian, while recomputing the 2327 nonlinear portion. 2328 2329 Collect on Mat 2330 2331 Input Parameters: 2332 . mat - the matrix (currently on AIJ matrices support this option) 2333 2334 Level: advanced 2335 2336 .seealso: MatStoreValues() 2337 2338 @*/ 2339 int MatRetrieveValues(Mat mat) 2340 { 2341 int ierr,(*f)(Mat); 2342 2343 PetscFunctionBegin; 2344 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2345 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2346 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2347 2348 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2349 if (f) { 2350 ierr = (*f)(mat);CHKERRQ(ierr); 2351 } else { 2352 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2353 } 2354 PetscFunctionReturn(0); 2355 } 2356 2357 /* 2358 This allows SeqAIJ matrices to be passed to the matlab engine 2359 */ 2360 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2361 #include "engine.h" /* Matlab include file */ 2362 #include "mex.h" /* Matlab include file */ 2363 EXTERN_C_BEGIN 2364 #undef __FUNCT__ 2365 #define __FUNCT__ "MatMatlabEnginePut_SeqAIJ" 2366 int MatMatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine) 2367 { 2368 int ierr,i,*ai,*aj; 2369 Mat B = (Mat)obj; 2370 mxArray *mat; 2371 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2372 2373 PetscFunctionBegin; 2374 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2375 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2376 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2377 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2378 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2379 2380 /* Matlab indices start at 0 for sparse (what a surprise) */ 2381 if (aij->indexshift) { 2382 ai = mxGetJc(mat); 2383 for (i=0; i<B->m+1; i++) { 2384 ai[i]--; 2385 } 2386 aj = mxGetIr(mat); 2387 for (i=0; i<aij->nz; i++) { 2388 aj[i]--; 2389 } 2390 } 2391 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2392 mxSetName(mat,obj->name); 2393 engPutArray((Engine *)mengine,mat); 2394 PetscFunctionReturn(0); 2395 } 2396 EXTERN_C_END 2397 2398 EXTERN_C_BEGIN 2399 #undef __FUNCT__ 2400 #define __FUNCT__ "MatMatlabEngineGet_SeqAIJ" 2401 int MatMatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine) 2402 { 2403 int ierr,ii; 2404 Mat mat = (Mat)obj; 2405 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2406 mxArray *mmat; 2407 2408 PetscFunctionBegin; 2409 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2410 aij->indexshift = 0; 2411 2412 mmat = engGetArray((Engine *)mengine,obj->name); 2413 2414 aij->nz = (mxGetJc(mmat))[mat->m]; 2415 ierr = PetscMalloc(((size_t) aij->nz)*(sizeof(int)+sizeof(PetscScalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 2416 aij->j = (int*)(aij->a + aij->nz); 2417 aij->i = aij->j + aij->nz; 2418 aij->singlemalloc = PETSC_TRUE; 2419 aij->freedata = PETSC_TRUE; 2420 2421 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2422 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2423 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 2424 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 2425 2426 for (ii=0; ii<mat->m; ii++) { 2427 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 2428 } 2429 2430 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2431 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2432 2433 PetscFunctionReturn(0); 2434 } 2435 EXTERN_C_END 2436 #endif 2437 2438 /* --------------------------------------------------------------------------------*/ 2439 #undef __FUNCT__ 2440 #define __FUNCT__ "MatCreateSeqAIJ" 2441 /*@C 2442 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2443 (the default parallel PETSc format). For good matrix assembly performance 2444 the user should preallocate the matrix storage by setting the parameter nz 2445 (or the array nnz). By setting these parameters accurately, performance 2446 during matrix assembly can be increased by more than a factor of 50. 2447 2448 Collective on MPI_Comm 2449 2450 Input Parameters: 2451 + comm - MPI communicator, set to PETSC_COMM_SELF 2452 . m - number of rows 2453 . n - number of columns 2454 . nz - number of nonzeros per row (same for all rows) 2455 - nnz - array containing the number of nonzeros in the various rows 2456 (possibly different for each row) or PETSC_NULL 2457 2458 Output Parameter: 2459 . A - the matrix 2460 2461 Notes: 2462 The AIJ format (also called the Yale sparse matrix format or 2463 compressed row storage), is fully compatible with standard Fortran 77 2464 storage. That is, the stored row and column indices can begin at 2465 either one (as in Fortran) or zero. See the users' manual for details. 2466 2467 Specify the preallocated storage with either nz or nnz (not both). 2468 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2469 allocation. For large problems you MUST preallocate memory or you 2470 will get TERRIBLE performance, see the users' manual chapter on matrices. 2471 2472 By default, this format uses inodes (identical nodes) when possible, to 2473 improve numerical efficiency of matrix-vector products and solves. We 2474 search for consecutive rows with the same nonzero structure, thereby 2475 reusing matrix information to achieve increased efficiency. 2476 2477 Options Database Keys: 2478 + -mat_aij_no_inode - Do not use inodes 2479 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2480 - -mat_aij_oneindex - Internally use indexing starting at 1 2481 rather than 0. Note that when calling MatSetValues(), 2482 the user still MUST index entries starting at 0! 2483 2484 Level: intermediate 2485 2486 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2487 2488 @*/ 2489 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz,Mat *A) 2490 { 2491 int ierr; 2492 2493 PetscFunctionBegin; 2494 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2495 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2496 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2497 PetscFunctionReturn(0); 2498 } 2499 2500 #define SKIP_ALLOCATION -4 2501 2502 #undef __FUNCT__ 2503 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2504 /*@C 2505 MatSeqAIJSetPreallocation - For good matrix assembly performance 2506 the user should preallocate the matrix storage by setting the parameter nz 2507 (or the array nnz). By setting these parameters accurately, performance 2508 during matrix assembly can be increased by more than a factor of 50. 2509 2510 Collective on MPI_Comm 2511 2512 Input Parameters: 2513 + comm - MPI communicator, set to PETSC_COMM_SELF 2514 . m - number of rows 2515 . n - number of columns 2516 . nz - number of nonzeros per row (same for all rows) 2517 - nnz - array containing the number of nonzeros in the various rows 2518 (possibly different for each row) or PETSC_NULL 2519 2520 Output Parameter: 2521 . A - the matrix 2522 2523 Notes: 2524 The AIJ format (also called the Yale sparse matrix format or 2525 compressed row storage), is fully compatible with standard Fortran 77 2526 storage. That is, the stored row and column indices can begin at 2527 either one (as in Fortran) or zero. See the users' manual for details. 2528 2529 Specify the preallocated storage with either nz or nnz (not both). 2530 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2531 allocation. For large problems you MUST preallocate memory or you 2532 will get TERRIBLE performance, see the users' manual chapter on matrices. 2533 2534 By default, this format uses inodes (identical nodes) when possible, to 2535 improve numerical efficiency of matrix-vector products and solves. We 2536 search for consecutive rows with the same nonzero structure, thereby 2537 reusing matrix information to achieve increased efficiency. 2538 2539 Options Database Keys: 2540 + -mat_aij_no_inode - Do not use inodes 2541 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2542 - -mat_aij_oneindex - Internally use indexing starting at 1 2543 rather than 0. Note that when calling MatSetValues(), 2544 the user still MUST index entries starting at 0! 2545 2546 Level: intermediate 2547 2548 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2549 2550 @*/ 2551 int MatSeqAIJSetPreallocation(Mat B,int nz,int *nnz) 2552 { 2553 Mat_SeqAIJ *b; 2554 size_t len = 0; 2555 PetscTruth flg2,skipallocation = PETSC_FALSE; 2556 int i,ierr; 2557 2558 PetscFunctionBegin; 2559 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flg2);CHKERRQ(ierr); 2560 if (!flg2) PetscFunctionReturn(0); 2561 2562 if (nz == SKIP_ALLOCATION) { 2563 skipallocation = PETSC_TRUE; 2564 nz = 0; 2565 } 2566 2567 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2568 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2569 if (nnz) { 2570 for (i=0; i<B->m; i++) { 2571 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2572 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); 2573 } 2574 } 2575 2576 B->preallocated = PETSC_TRUE; 2577 b = (Mat_SeqAIJ*)B->data; 2578 2579 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2580 if (!nnz) { 2581 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2582 else if (nz <= 0) nz = 1; 2583 for (i=0; i<B->m; i++) b->imax[i] = nz; 2584 nz = nz*B->m; 2585 } else { 2586 nz = 0; 2587 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2588 } 2589 2590 if (!skipallocation) { 2591 /* allocate the matrix space */ 2592 len = ((size_t) nz)*(sizeof(int) + sizeof(PetscScalar)) + (B->m+1)*sizeof(int); 2593 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2594 b->j = (int*)(b->a + nz); 2595 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2596 b->i = b->j + nz; 2597 b->i[0] = -b->indexshift; 2598 for (i=1; i<B->m+1; i++) { 2599 b->i[i] = b->i[i-1] + b->imax[i-1]; 2600 } 2601 b->singlemalloc = PETSC_TRUE; 2602 b->freedata = PETSC_TRUE; 2603 } else { 2604 b->freedata = PETSC_FALSE; 2605 } 2606 2607 /* b->ilen will count nonzeros in each row so far. */ 2608 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2609 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2610 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2611 2612 b->nz = 0; 2613 b->maxnz = nz; 2614 B->info.nz_unneeded = (double)b->maxnz; 2615 PetscFunctionReturn(0); 2616 } 2617 2618 EXTERN int RegisterApplyPtAPRoutines_Private(Mat); 2619 2620 EXTERN_C_BEGIN 2621 extern int MatConvert_SeqAIJ_SeqSBAIJ(Mat,MatType,Mat*); 2622 extern int MatConvert_SeqAIJ_SeqBAIJ(Mat,MatType,Mat*); 2623 EXTERN_C_END 2624 2625 EXTERN_C_BEGIN 2626 #undef __FUNCT__ 2627 #define __FUNCT__ "MatCreate_SeqAIJ" 2628 int MatCreate_SeqAIJ(Mat B) 2629 { 2630 Mat_SeqAIJ *b; 2631 int ierr,size; 2632 PetscTruth flg; 2633 2634 PetscFunctionBegin; 2635 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2636 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2637 2638 B->m = B->M = PetscMax(B->m,B->M); 2639 B->n = B->N = PetscMax(B->n,B->N); 2640 2641 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2642 B->data = (void*)b; 2643 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2644 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2645 B->factor = 0; 2646 B->lupivotthreshold = 1.0; 2647 B->mapping = 0; 2648 ierr = PetscOptionsGetReal(B->prefix,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2649 ierr = PetscOptionsHasName(B->prefix,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2650 b->row = 0; 2651 b->col = 0; 2652 b->icol = 0; 2653 b->indexshift = 0; 2654 b->reallocs = 0; 2655 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_oneindex",&flg);CHKERRQ(ierr); 2656 if (flg) b->indexshift = -1; 2657 2658 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2659 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2660 2661 b->sorted = PETSC_FALSE; 2662 b->ignorezeroentries = PETSC_FALSE; 2663 b->roworiented = PETSC_TRUE; 2664 b->nonew = 0; 2665 b->diag = 0; 2666 b->solve_work = 0; 2667 B->spptr = 0; 2668 b->inode.use = PETSC_TRUE; 2669 b->inode.node_count = 0; 2670 b->inode.size = 0; 2671 b->inode.limit = 5; 2672 b->inode.max_limit = 5; 2673 b->saved_values = 0; 2674 b->idiag = 0; 2675 b->ssor = 0; 2676 b->keepzeroedrows = PETSC_FALSE; 2677 2678 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2679 2680 #if defined(PETSC_HAVE_SUPERLU) 2681 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_superlu",&flg);CHKERRQ(ierr); 2682 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B);CHKERRQ(ierr); } 2683 #endif 2684 2685 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_essl",&flg);CHKERRQ(ierr); 2686 if (flg) { ierr = MatUseEssl_SeqAIJ(B);CHKERRQ(ierr); } 2687 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_lusol",&flg);CHKERRQ(ierr); 2688 if (flg) { ierr = MatUseLUSOL_SeqAIJ(B);CHKERRQ(ierr); } 2689 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_matlab",&flg);CHKERRQ(ierr); 2690 if (flg) {ierr = MatUseMatlab_SeqAIJ(B);CHKERRQ(ierr);} 2691 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_dxml",&flg);CHKERRQ(ierr); 2692 if (flg) { 2693 if (!b->indexshift) SETERRQ(PETSC_ERR_LIB,"need -mat_aij_oneindex with -mat_aij_dxml"); 2694 ierr = MatUseDXML_SeqAIJ(B);CHKERRQ(ierr); 2695 } 2696 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2697 "MatSeqAIJSetColumnIndices_SeqAIJ", 2698 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2699 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2700 "MatStoreValues_SeqAIJ", 2701 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2702 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2703 "MatRetrieveValues_SeqAIJ", 2704 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2705 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2706 "MatConvert_SeqAIJ_SeqSBAIJ", 2707 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2708 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2709 "MatConvert_SeqAIJ_SeqBAIJ", 2710 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2711 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2712 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2713 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2714 #endif 2715 ierr = RegisterApplyPtAPRoutines_Private(B);CHKERRQ(ierr); 2716 PetscFunctionReturn(0); 2717 } 2718 EXTERN_C_END 2719 2720 #undef __FUNCT__ 2721 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2722 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2723 { 2724 Mat C; 2725 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2726 int i,m = A->m,shift = a->indexshift,ierr; 2727 size_t len; 2728 2729 PetscFunctionBegin; 2730 *B = 0; 2731 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2732 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 2733 c = (Mat_SeqAIJ*)C->data; 2734 2735 C->factor = A->factor; 2736 c->row = 0; 2737 c->col = 0; 2738 c->icol = 0; 2739 c->indexshift = shift; 2740 c->keepzeroedrows = a->keepzeroedrows; 2741 C->assembled = PETSC_TRUE; 2742 2743 C->M = A->m; 2744 C->N = A->n; 2745 2746 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2747 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2748 for (i=0; i<m; i++) { 2749 c->imax[i] = a->imax[i]; 2750 c->ilen[i] = a->ilen[i]; 2751 } 2752 2753 /* allocate the matrix space */ 2754 c->singlemalloc = PETSC_TRUE; 2755 len = ((size_t) (m+1))*sizeof(int)+(a->i[m])*(sizeof(PetscScalar)+sizeof(int)); 2756 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2757 c->j = (int*)(c->a + a->i[m] + shift); 2758 c->i = c->j + a->i[m] + shift; 2759 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2760 if (m > 0) { 2761 ierr = PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int));CHKERRQ(ierr); 2762 if (cpvalues == MAT_COPY_VALUES) { 2763 ierr = PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(PetscScalar));CHKERRQ(ierr); 2764 } else { 2765 ierr = PetscMemzero(c->a,(a->i[m]+shift)*sizeof(PetscScalar));CHKERRQ(ierr); 2766 } 2767 } 2768 2769 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2770 c->sorted = a->sorted; 2771 c->roworiented = a->roworiented; 2772 c->nonew = a->nonew; 2773 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2774 c->saved_values = 0; 2775 c->idiag = 0; 2776 c->ssor = 0; 2777 c->ignorezeroentries = a->ignorezeroentries; 2778 c->freedata = PETSC_TRUE; 2779 2780 if (a->diag) { 2781 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2782 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2783 for (i=0; i<m; i++) { 2784 c->diag[i] = a->diag[i]; 2785 } 2786 } else c->diag = 0; 2787 c->inode.use = a->inode.use; 2788 c->inode.limit = a->inode.limit; 2789 c->inode.max_limit = a->inode.max_limit; 2790 if (a->inode.size){ 2791 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2792 c->inode.node_count = a->inode.node_count; 2793 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2794 } else { 2795 c->inode.size = 0; 2796 c->inode.node_count = 0; 2797 } 2798 c->nz = a->nz; 2799 c->maxnz = a->maxnz; 2800 c->solve_work = 0; 2801 C->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2802 C->preallocated = PETSC_TRUE; 2803 2804 *B = C; 2805 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2806 PetscFunctionReturn(0); 2807 } 2808 2809 EXTERN_C_BEGIN 2810 #undef __FUNCT__ 2811 #define __FUNCT__ "MatLoad_SeqAIJ" 2812 int MatLoad_SeqAIJ(PetscViewer viewer,MatType type,Mat *A) 2813 { 2814 Mat_SeqAIJ *a; 2815 Mat B; 2816 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N,shift; 2817 MPI_Comm comm; 2818 #if defined(PETSC_HAVE_SPOOLES) || defined(PETSC_HAVE_SUPERLU) || defined(PETSC_HAVE_SUPERLUDIST) || defined(PETSC_HAVE_UMFPACK) 2819 PetscTruth flag; 2820 #endif 2821 2822 PetscFunctionBegin; 2823 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2824 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2825 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2826 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2827 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2828 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2829 M = header[1]; N = header[2]; nz = header[3]; 2830 2831 if (nz < 0) { 2832 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2833 } 2834 2835 /* read in row lengths */ 2836 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2837 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2838 2839 /* create our matrix */ 2840 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A);CHKERRQ(ierr); 2841 B = *A; 2842 a = (Mat_SeqAIJ*)B->data; 2843 shift = a->indexshift; 2844 2845 /* read in column indices and adjust for Fortran indexing*/ 2846 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2847 if (shift) { 2848 for (i=0; i<nz; i++) { 2849 a->j[i] += 1; 2850 } 2851 } 2852 2853 /* read in nonzero values */ 2854 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2855 2856 /* set matrix "i" values */ 2857 a->i[0] = -shift; 2858 for (i=1; i<= M; i++) { 2859 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2860 a->ilen[i-1] = rowlengths[i-1]; 2861 } 2862 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2863 2864 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2865 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2866 #if defined(PETSC_HAVE_SPOOLES) 2867 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_spooles",&flag);CHKERRQ(ierr); 2868 if (flag) { ierr = MatUseSpooles_SeqAIJ(B);CHKERRQ(ierr); } 2869 #endif 2870 #if defined(PETSC_HAVE_SUPERLU) 2871 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_superlu",&flag);CHKERRQ(ierr); 2872 if (flag) { ierr = MatUseSuperLU_SeqAIJ(B);CHKERRQ(ierr); } 2873 #endif 2874 #if defined(PETSC_HAVE_SUPERLUDIST) 2875 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_superlu_dist",&flag);CHKERRQ(ierr); 2876 if (flag) { ierr = MatUseSuperLU_DIST_MPIAIJ(B);CHKERRQ(ierr); } 2877 #endif 2878 #if defined(PETSC_HAVE_UMFPACK) 2879 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_umfpack",&flag);CHKERRQ(ierr); 2880 if (flag) { ierr = MatUseUMFPACK_SeqAIJ(B);CHKERRQ(ierr); } 2881 #endif 2882 PetscFunctionReturn(0); 2883 } 2884 EXTERN_C_END 2885 2886 #undef __FUNCT__ 2887 #define __FUNCT__ "MatEqual_SeqAIJ" 2888 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2889 { 2890 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2891 int ierr; 2892 PetscTruth flag; 2893 2894 PetscFunctionBegin; 2895 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flag);CHKERRQ(ierr); 2896 if (!flag) SETERRQ(PETSC_ERR_ARG_INCOMP,"Matrices must be same type"); 2897 2898 /* If the matrix dimensions are not equal,or no of nonzeros or shift */ 2899 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)|| (a->indexshift != b->indexshift)) { 2900 *flg = PETSC_FALSE; 2901 PetscFunctionReturn(0); 2902 } 2903 2904 /* if the a->i are the same */ 2905 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2906 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2907 2908 /* if a->j are the same */ 2909 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2910 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2911 2912 /* if a->a are the same */ 2913 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2914 2915 PetscFunctionReturn(0); 2916 2917 } 2918 2919 #undef __FUNCT__ 2920 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2921 /*@C 2922 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2923 provided by the user. 2924 2925 Coolective on MPI_Comm 2926 2927 Input Parameters: 2928 + comm - must be an MPI communicator of size 1 2929 . m - number of rows 2930 . n - number of columns 2931 . i - row indices 2932 . j - column indices 2933 - a - matrix values 2934 2935 Output Parameter: 2936 . mat - the matrix 2937 2938 Level: intermediate 2939 2940 Notes: 2941 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2942 once the matrix is destroyed 2943 2944 You cannot set new nonzero locations into this matrix, that will generate an error. 2945 2946 The i and j indices can be either 0- or 1 based 2947 2948 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2949 2950 @*/ 2951 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,PetscScalar *a,Mat *mat) 2952 { 2953 int ierr,ii; 2954 Mat_SeqAIJ *aij; 2955 2956 PetscFunctionBegin; 2957 ierr = MatCreateSeqAIJ(comm,m,n,SKIP_ALLOCATION,0,mat);CHKERRQ(ierr); 2958 aij = (Mat_SeqAIJ*)(*mat)->data; 2959 2960 if (i[0] == 1) { 2961 aij->indexshift = -1; 2962 } else if (i[0]) { 2963 SETERRQ(1,"i (row indices) do not start with 0 or 1"); 2964 } 2965 aij->i = i; 2966 aij->j = j; 2967 aij->a = a; 2968 aij->singlemalloc = PETSC_FALSE; 2969 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2970 aij->freedata = PETSC_FALSE; 2971 2972 for (ii=0; ii<m; ii++) { 2973 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2974 #if defined(PETSC_USE_BOPT_g) 2975 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]); 2976 #endif 2977 } 2978 #if defined(PETSC_USE_BOPT_g) 2979 for (ii=0; ii<aij->i[m]; ii++) { 2980 if (j[ii] < -aij->indexshift) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 2981 if (j[ii] > n - 1 -aij->indexshift) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 2982 } 2983 #endif 2984 2985 /* changes indices to start at 0 */ 2986 if (i[0]) { 2987 aij->indexshift = 0; 2988 for (ii=0; ii<m; ii++) { 2989 i[ii]--; 2990 } 2991 for (ii=0; ii<i[m]; ii++) { 2992 j[ii]--; 2993 } 2994 } 2995 2996 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2997 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2998 PetscFunctionReturn(0); 2999 } 3000 3001 #undef __FUNCT__ 3002 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3003 int MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3004 { 3005 int ierr; 3006 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3007 3008 PetscFunctionBegin; 3009 if (coloring->ctype == IS_COLORING_LOCAL) { 3010 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3011 a->coloring = coloring; 3012 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3013 int i,*larray; 3014 ISColoring ocoloring; 3015 ISColoringValue *colors; 3016 3017 /* set coloring for diagonal portion */ 3018 ierr = PetscMalloc((A->n+1)*sizeof(int),&larray);CHKERRQ(ierr); 3019 for (i=0; i<A->n; i++) { 3020 larray[i] = i; 3021 } 3022 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3023 ierr = PetscMalloc((A->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3024 for (i=0; i<A->n; i++) { 3025 colors[i] = coloring->colors[larray[i]]; 3026 } 3027 ierr = PetscFree(larray);CHKERRQ(ierr); 3028 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 3029 a->coloring = ocoloring; 3030 } 3031 PetscFunctionReturn(0); 3032 } 3033 3034 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 3035 EXTERN_C_BEGIN 3036 #include "adic/ad_utils.h" 3037 EXTERN_C_END 3038 3039 #undef __FUNCT__ 3040 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3041 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3042 { 3043 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3044 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3045 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3046 ISColoringValue *color; 3047 3048 PetscFunctionBegin; 3049 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3050 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3051 color = a->coloring->colors; 3052 /* loop over rows */ 3053 for (i=0; i<m; i++) { 3054 nz = ii[i+1] - ii[i]; 3055 /* loop over columns putting computed value into matrix */ 3056 for (j=0; j<nz; j++) { 3057 *v++ = values[color[*jj++]]; 3058 } 3059 values += nlen; /* jump to next row of derivatives */ 3060 } 3061 PetscFunctionReturn(0); 3062 } 3063 3064 #else 3065 3066 #undef __FUNCT__ 3067 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3068 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3069 { 3070 PetscFunctionBegin; 3071 SETERRQ(1,"PETSc installed without ADIC"); 3072 } 3073 3074 #endif 3075 3076 #undef __FUNCT__ 3077 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3078 int MatSetValuesAdifor_SeqAIJ(Mat A,int nl,void *advalues) 3079 { 3080 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3081 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j; 3082 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3083 ISColoringValue *color; 3084 3085 PetscFunctionBegin; 3086 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3087 color = a->coloring->colors; 3088 /* loop over rows */ 3089 for (i=0; i<m; i++) { 3090 nz = ii[i+1] - ii[i]; 3091 /* loop over columns putting computed value into matrix */ 3092 for (j=0; j<nz; j++) { 3093 *v++ = values[color[*jj++]]; 3094 } 3095 values += nl; /* jump to next row of derivatives */ 3096 } 3097 PetscFunctionReturn(0); 3098 } 3099 3100