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