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