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