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