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