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