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