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 int *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1446 int ma,na,mb,nb, i,ierr; 1447 1448 PetscFunctionBegin; 1449 bij = (Mat_SeqAIJ *) B->data; 1450 1451 ierr = MatGetSize(A,&ma,&na); CHKERRQ(ierr); 1452 ierr = MatGetSize(B,&mb,&nb); CHKERRQ(ierr); 1453 if (ma!=nb || na!=mb) 1454 SETERRQ(1,"Incompatible A/B sizes for symmetry test"); 1455 aii = aij->i; bii = bij->i; 1456 adx = aij->j; bdx = bij->j; 1457 va = aij->a; vb = bij->a; 1458 ierr = PetscMalloc(ma*sizeof(int),&aptr); CHKERRQ(ierr); 1459 ierr = PetscMalloc(mb*sizeof(int),&bptr); CHKERRQ(ierr); 1460 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1461 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1462 1463 *f = PETSC_TRUE; 1464 for (i=0; i<ma; i++) { 1465 /*printf("row %d spans %d--%d; we start @ %d\n", 1466 i,idx[ii[i]],idx[ii[i+1]-1],idx[aptr[i]]);*/ 1467 while (aptr[i]<aii[i+1]) { 1468 int idc,idr; PetscScalar vc,vr; 1469 /* column/row index/value */ 1470 idc = adx[aptr[i]]; idr = bdx[bptr[idc]]; 1471 vc = va[aptr[i]]; vr = vb[bptr[idc]]; 1472 /*printf("comparing %d: (%d,%d)=%e to (%d,%d)=%e\n", 1473 aptr[i],i,idc,vc,idc,idr,vr);*/ 1474 if (i!=idr || vc!=vr) { 1475 *f = PETSC_FALSE; goto done; 1476 } else { 1477 aptr[i]++; if (B || i!=idc) bptr[idc]++; 1478 } 1479 } 1480 } 1481 done: 1482 ierr = PetscFree(aptr); CHKERRQ(ierr); 1483 if (B) { 1484 ierr = PetscFree(bptr); CHKERRQ(ierr); 1485 } 1486 1487 PetscFunctionReturn(0); 1488 } 1489 EXTERN_C_END 1490 1491 #undef __FUNCT__ 1492 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1493 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1494 { 1495 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1496 PetscScalar *l,*r,x,*v; 1497 int ierr,i,j,m = A->m,n = A->n,M,nz = a->nz,*jj; 1498 1499 PetscFunctionBegin; 1500 if (ll) { 1501 /* The local size is used so that VecMPI can be passed to this routine 1502 by MatDiagonalScale_MPIAIJ */ 1503 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1504 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1505 ierr = VecGetArrayFast(ll,&l);CHKERRQ(ierr); 1506 v = a->a; 1507 for (i=0; i<m; i++) { 1508 x = l[i]; 1509 M = a->i[i+1] - a->i[i]; 1510 for (j=0; j<M; j++) { (*v++) *= x;} 1511 } 1512 ierr = VecRestoreArrayFast(ll,&l);CHKERRQ(ierr); 1513 PetscLogFlops(nz); 1514 } 1515 if (rr) { 1516 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1517 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1518 ierr = VecGetArrayFast(rr,&r);CHKERRQ(ierr); 1519 v = a->a; jj = a->j; 1520 for (i=0; i<nz; i++) { 1521 (*v++) *= r[*jj++]; 1522 } 1523 ierr = VecRestoreArrayFast(rr,&r);CHKERRQ(ierr); 1524 PetscLogFlops(nz); 1525 } 1526 PetscFunctionReturn(0); 1527 } 1528 1529 #undef __FUNCT__ 1530 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1531 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatReuse scall,Mat *B) 1532 { 1533 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1534 int *smap,i,k,kstart,kend,ierr,oldcols = A->n,*lens; 1535 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1536 int *irow,*icol,nrows,ncols; 1537 int *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1538 PetscScalar *a_new,*mat_a; 1539 Mat C; 1540 PetscTruth stride; 1541 1542 PetscFunctionBegin; 1543 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1544 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1545 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1546 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1547 1548 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1549 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1550 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1551 1552 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1553 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1554 if (stride && step == 1) { 1555 /* special case of contiguous rows */ 1556 ierr = PetscMalloc((2*nrows+1)*sizeof(int),&lens);CHKERRQ(ierr); 1557 starts = lens + nrows; 1558 /* loop over new rows determining lens and starting points */ 1559 for (i=0; i<nrows; i++) { 1560 kstart = ai[irow[i]]; 1561 kend = kstart + ailen[irow[i]]; 1562 for (k=kstart; k<kend; k++) { 1563 if (aj[k] >= first) { 1564 starts[i] = k; 1565 break; 1566 } 1567 } 1568 sum = 0; 1569 while (k < kend) { 1570 if (aj[k++] >= first+ncols) break; 1571 sum++; 1572 } 1573 lens[i] = sum; 1574 } 1575 /* create submatrix */ 1576 if (scall == MAT_REUSE_MATRIX) { 1577 int n_cols,n_rows; 1578 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1579 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1580 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1581 C = *B; 1582 } else { 1583 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1584 } 1585 c = (Mat_SeqAIJ*)C->data; 1586 1587 /* loop over rows inserting into submatrix */ 1588 a_new = c->a; 1589 j_new = c->j; 1590 i_new = c->i; 1591 1592 for (i=0; i<nrows; i++) { 1593 ii = starts[i]; 1594 lensi = lens[i]; 1595 for (k=0; k<lensi; k++) { 1596 *j_new++ = aj[ii+k] - first; 1597 } 1598 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1599 a_new += lensi; 1600 i_new[i+1] = i_new[i] + lensi; 1601 c->ilen[i] = lensi; 1602 } 1603 ierr = PetscFree(lens);CHKERRQ(ierr); 1604 } else { 1605 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1606 ierr = PetscMalloc((1+oldcols)*sizeof(int),&smap);CHKERRQ(ierr); 1607 1608 ierr = PetscMalloc((1+nrows)*sizeof(int),&lens);CHKERRQ(ierr); 1609 ierr = PetscMemzero(smap,oldcols*sizeof(int));CHKERRQ(ierr); 1610 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1611 /* determine lens of each row */ 1612 for (i=0; i<nrows; i++) { 1613 kstart = ai[irow[i]]; 1614 kend = kstart + a->ilen[irow[i]]; 1615 lens[i] = 0; 1616 for (k=kstart; k<kend; k++) { 1617 if (smap[aj[k]]) { 1618 lens[i]++; 1619 } 1620 } 1621 } 1622 /* Create and fill new matrix */ 1623 if (scall == MAT_REUSE_MATRIX) { 1624 PetscTruth equal; 1625 1626 c = (Mat_SeqAIJ *)((*B)->data); 1627 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1628 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(int),&equal);CHKERRQ(ierr); 1629 if (!equal) { 1630 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1631 } 1632 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(int));CHKERRQ(ierr); 1633 C = *B; 1634 } else { 1635 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1636 } 1637 c = (Mat_SeqAIJ *)(C->data); 1638 for (i=0; i<nrows; i++) { 1639 row = irow[i]; 1640 kstart = ai[row]; 1641 kend = kstart + a->ilen[row]; 1642 mat_i = c->i[i]; 1643 mat_j = c->j + mat_i; 1644 mat_a = c->a + mat_i; 1645 mat_ilen = c->ilen + i; 1646 for (k=kstart; k<kend; k++) { 1647 if ((tcol=smap[a->j[k]])) { 1648 *mat_j++ = tcol - 1; 1649 *mat_a++ = a->a[k]; 1650 (*mat_ilen)++; 1651 1652 } 1653 } 1654 } 1655 /* Free work space */ 1656 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1657 ierr = PetscFree(smap);CHKERRQ(ierr); 1658 ierr = PetscFree(lens);CHKERRQ(ierr); 1659 } 1660 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1661 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1662 1663 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1664 *B = C; 1665 PetscFunctionReturn(0); 1666 } 1667 1668 /* 1669 */ 1670 #undef __FUNCT__ 1671 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1672 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1673 { 1674 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1675 int ierr; 1676 Mat outA; 1677 PetscTruth row_identity,col_identity; 1678 1679 PetscFunctionBegin; 1680 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1681 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1682 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1683 if (!row_identity || !col_identity) { 1684 SETERRQ(1,"Row and column permutations must be identity for in-place ILU"); 1685 } 1686 1687 outA = inA; 1688 inA->factor = FACTOR_LU; 1689 a->row = row; 1690 a->col = col; 1691 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1692 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1693 1694 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1695 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1696 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1697 PetscLogObjectParent(inA,a->icol); 1698 1699 if (!a->solve_work) { /* this matrix may have been factored before */ 1700 ierr = PetscMalloc((inA->m+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1701 } 1702 1703 if (!a->diag) { 1704 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1705 } 1706 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1707 PetscFunctionReturn(0); 1708 } 1709 1710 #include "petscblaslapack.h" 1711 #undef __FUNCT__ 1712 #define __FUNCT__ "MatScale_SeqAIJ" 1713 int MatScale_SeqAIJ(PetscScalar *alpha,Mat inA) 1714 { 1715 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1716 int one = 1; 1717 1718 PetscFunctionBegin; 1719 BLscal_(&a->nz,alpha,a->a,&one); 1720 PetscLogFlops(a->nz); 1721 PetscFunctionReturn(0); 1722 } 1723 1724 #undef __FUNCT__ 1725 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1726 int MatGetSubMatrices_SeqAIJ(Mat A,int n,IS *irow,IS *icol,MatReuse scall,Mat **B) 1727 { 1728 int ierr,i; 1729 1730 PetscFunctionBegin; 1731 if (scall == MAT_INITIAL_MATRIX) { 1732 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1733 } 1734 1735 for (i=0; i<n; i++) { 1736 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1737 } 1738 PetscFunctionReturn(0); 1739 } 1740 1741 #undef __FUNCT__ 1742 #define __FUNCT__ "MatGetBlockSize_SeqAIJ" 1743 int MatGetBlockSize_SeqAIJ(Mat A,int *bs) 1744 { 1745 PetscFunctionBegin; 1746 *bs = 1; 1747 PetscFunctionReturn(0); 1748 } 1749 1750 #undef __FUNCT__ 1751 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1752 int MatIncreaseOverlap_SeqAIJ(Mat A,int is_max,IS *is,int ov) 1753 { 1754 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1755 int row,i,j,k,l,m,n,*idx,ierr,*nidx,isz,val; 1756 int start,end,*ai,*aj; 1757 PetscBT table; 1758 1759 PetscFunctionBegin; 1760 m = A->m; 1761 ai = a->i; 1762 aj = a->j; 1763 1764 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal overlap value used"); 1765 1766 ierr = PetscMalloc((m+1)*sizeof(int),&nidx);CHKERRQ(ierr); 1767 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1768 1769 for (i=0; i<is_max; i++) { 1770 /* Initialize the two local arrays */ 1771 isz = 0; 1772 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1773 1774 /* Extract the indices, assume there can be duplicate entries */ 1775 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1776 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1777 1778 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1779 for (j=0; j<n ; ++j){ 1780 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1781 } 1782 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1783 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1784 1785 k = 0; 1786 for (j=0; j<ov; j++){ /* for each overlap */ 1787 n = isz; 1788 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1789 row = nidx[k]; 1790 start = ai[row]; 1791 end = ai[row+1]; 1792 for (l = start; l<end ; l++){ 1793 val = aj[l] ; 1794 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1795 } 1796 } 1797 } 1798 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1799 } 1800 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1801 ierr = PetscFree(nidx);CHKERRQ(ierr); 1802 PetscFunctionReturn(0); 1803 } 1804 1805 /* -------------------------------------------------------------- */ 1806 #undef __FUNCT__ 1807 #define __FUNCT__ "MatPermute_SeqAIJ" 1808 int MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1809 { 1810 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1811 PetscScalar *vwork; 1812 int i,ierr,nz,m = A->m,n = A->n,*cwork; 1813 int *row,*col,*cnew,j,*lens; 1814 IS icolp,irowp; 1815 1816 PetscFunctionBegin; 1817 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1818 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1819 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1820 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1821 1822 /* determine lengths of permuted rows */ 1823 ierr = PetscMalloc((m+1)*sizeof(int),&lens);CHKERRQ(ierr); 1824 for (i=0; i<m; i++) { 1825 lens[row[i]] = a->i[i+1] - a->i[i]; 1826 } 1827 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1828 ierr = PetscFree(lens);CHKERRQ(ierr); 1829 1830 ierr = PetscMalloc(n*sizeof(int),&cnew);CHKERRQ(ierr); 1831 for (i=0; i<m; i++) { 1832 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1833 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1834 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1835 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1836 } 1837 ierr = PetscFree(cnew);CHKERRQ(ierr); 1838 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1839 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1840 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1841 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1842 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1843 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1844 PetscFunctionReturn(0); 1845 } 1846 1847 #undef __FUNCT__ 1848 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1849 int MatPrintHelp_SeqAIJ(Mat A) 1850 { 1851 static PetscTruth called = PETSC_FALSE; 1852 MPI_Comm comm = A->comm; 1853 int ierr; 1854 1855 PetscFunctionBegin; 1856 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1857 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1858 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1859 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1860 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1861 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1862 #if defined(PETSC_HAVE_ESSL) 1863 ierr = (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n");CHKERRQ(ierr); 1864 #endif 1865 #if defined(PETSC_HAVE_LUSOL) 1866 ierr = (*PetscHelpPrintf)(comm," -mat_aij_lusol: Use the Stanford LUSOL sparse factorization and solve.\n");CHKERRQ(ierr); 1867 #endif 1868 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 1869 ierr = (*PetscHelpPrintf)(comm," -mat_aij_matlab: Use Matlab engine sparse LU factorization and solve.\n");CHKERRQ(ierr); 1870 #endif 1871 PetscFunctionReturn(0); 1872 } 1873 EXTERN int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg); 1874 EXTERN int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1875 EXTERN int MatILUDTFactor_SeqAIJ(Mat,MatFactorInfo*,IS,IS,Mat*); 1876 #undef __FUNCT__ 1877 #define __FUNCT__ "MatCopy_SeqAIJ" 1878 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1879 { 1880 int ierr; 1881 PetscTruth flg; 1882 1883 PetscFunctionBegin; 1884 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flg);CHKERRQ(ierr); 1885 if (str == SAME_NONZERO_PATTERN && flg) { 1886 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1887 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1888 1889 if (a->i[A->m] != b->i[B->m]) { 1890 SETERRQ(1,"Number of nonzeros in two matrices are different"); 1891 } 1892 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 1893 } else { 1894 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1895 } 1896 PetscFunctionReturn(0); 1897 } 1898 1899 #undef __FUNCT__ 1900 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1901 int MatSetUpPreallocation_SeqAIJ(Mat A) 1902 { 1903 int ierr; 1904 1905 PetscFunctionBegin; 1906 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1907 PetscFunctionReturn(0); 1908 } 1909 1910 #undef __FUNCT__ 1911 #define __FUNCT__ "MatGetArray_SeqAIJ" 1912 int MatGetArray_SeqAIJ(Mat A,PetscScalar **array) 1913 { 1914 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1915 PetscFunctionBegin; 1916 *array = a->a; 1917 PetscFunctionReturn(0); 1918 } 1919 1920 #undef __FUNCT__ 1921 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1922 int MatRestoreArray_SeqAIJ(Mat A,PetscScalar **array) 1923 { 1924 PetscFunctionBegin; 1925 PetscFunctionReturn(0); 1926 } 1927 1928 #undef __FUNCT__ 1929 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1930 int MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1931 { 1932 int (*f)(void *,Vec,Vec,void*) = (int (*)(void *,Vec,Vec,void *))coloring->f; 1933 int k,ierr,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1934 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1935 PetscScalar *vscale_array; 1936 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1937 Vec w1,w2,w3; 1938 void *fctx = coloring->fctx; 1939 PetscTruth flg; 1940 1941 PetscFunctionBegin; 1942 if (!coloring->w1) { 1943 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1944 PetscLogObjectParent(coloring,coloring->w1); 1945 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1946 PetscLogObjectParent(coloring,coloring->w2); 1947 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1948 PetscLogObjectParent(coloring,coloring->w3); 1949 } 1950 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1951 1952 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1953 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1954 if (flg) { 1955 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 1956 } else { 1957 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1958 } 1959 1960 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1961 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1962 1963 /* 1964 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1965 coloring->F for the coarser grids from the finest 1966 */ 1967 if (coloring->F) { 1968 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1969 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1970 if (m1 != m2) { 1971 coloring->F = 0; 1972 } 1973 } 1974 1975 if (coloring->F) { 1976 w1 = coloring->F; 1977 coloring->F = 0; 1978 } else { 1979 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1980 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1981 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1982 } 1983 1984 /* 1985 Compute all the scale factors and share with other processors 1986 */ 1987 ierr = VecGetArrayFast(x1,&xx);CHKERRQ(ierr);xx = xx - start; 1988 ierr = VecGetArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 1989 for (k=0; k<coloring->ncolors; k++) { 1990 /* 1991 Loop over each column associated with color adding the 1992 perturbation to the vector w3. 1993 */ 1994 for (l=0; l<coloring->ncolumns[k]; l++) { 1995 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1996 dx = xx[col]; 1997 if (dx == 0.0) dx = 1.0; 1998 #if !defined(PETSC_USE_COMPLEX) 1999 if (dx < umin && dx >= 0.0) dx = umin; 2000 else if (dx < 0.0 && dx > -umin) dx = -umin; 2001 #else 2002 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2003 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2004 #endif 2005 dx *= epsilon; 2006 vscale_array[col] = 1.0/dx; 2007 } 2008 } 2009 vscale_array = vscale_array + start;ierr = VecRestoreArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2010 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2011 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2012 2013 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2014 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2015 2016 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2017 else vscaleforrow = coloring->columnsforrow; 2018 2019 ierr = VecGetArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2020 /* 2021 Loop over each color 2022 */ 2023 for (k=0; k<coloring->ncolors; k++) { 2024 coloring->currentcolor = k; 2025 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2026 ierr = VecGetArrayFast(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2027 /* 2028 Loop over each column associated with color adding the 2029 perturbation to the vector w3. 2030 */ 2031 for (l=0; l<coloring->ncolumns[k]; l++) { 2032 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2033 dx = xx[col]; 2034 if (dx == 0.0) dx = 1.0; 2035 #if !defined(PETSC_USE_COMPLEX) 2036 if (dx < umin && dx >= 0.0) dx = umin; 2037 else if (dx < 0.0 && dx > -umin) dx = -umin; 2038 #else 2039 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2040 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2041 #endif 2042 dx *= epsilon; 2043 if (!PetscAbsScalar(dx)) SETERRQ(1,"Computed 0 differencing parameter"); 2044 w3_array[col] += dx; 2045 } 2046 w3_array = w3_array + start; ierr = VecRestoreArrayFast(w3,&w3_array);CHKERRQ(ierr); 2047 2048 /* 2049 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2050 */ 2051 2052 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2053 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2054 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2055 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 2056 2057 /* 2058 Loop over rows of vector, putting results into Jacobian matrix 2059 */ 2060 ierr = VecGetArrayFast(w2,&y);CHKERRQ(ierr); 2061 for (l=0; l<coloring->nrows[k]; l++) { 2062 row = coloring->rows[k][l]; 2063 col = coloring->columnsforrow[k][l]; 2064 y[row] *= vscale_array[vscaleforrow[k][l]]; 2065 srow = row + start; 2066 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2067 } 2068 ierr = VecRestoreArrayFast(w2,&y);CHKERRQ(ierr); 2069 } 2070 coloring->currentcolor = k; 2071 ierr = VecRestoreArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2072 xx = xx + start; ierr = VecRestoreArrayFast(x1,&xx);CHKERRQ(ierr); 2073 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2074 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2075 PetscFunctionReturn(0); 2076 } 2077 2078 #include "petscblaslapack.h" 2079 #undef __FUNCT__ 2080 #define __FUNCT__ "MatAXPY_SeqAIJ" 2081 int MatAXPY_SeqAIJ(PetscScalar *a,Mat X,Mat Y,MatStructure str) 2082 { 2083 int ierr,one=1,i; 2084 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2085 2086 PetscFunctionBegin; 2087 if (str == SAME_NONZERO_PATTERN) { 2088 BLaxpy_(&x->nz,a,x->a,&one,y->a,&one); 2089 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2090 if (y->xtoy && y->XtoY != X) { 2091 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2092 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2093 } 2094 if (!y->xtoy) { /* get xtoy */ 2095 ierr = MatAXPYGetxtoy_Private(X->m,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2096 y->XtoY = X; 2097 } 2098 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += (*a)*(x->a[i]); 2099 PetscLogInfo(0,"MatAXPY_SeqAIJ: ratio of nnz(X)/nnz(Y): %d/%d = %g\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz); 2100 } else { 2101 ierr = MatAXPY_Basic(a,X,Y,str);CHKERRQ(ierr); 2102 } 2103 PetscFunctionReturn(0); 2104 } 2105 2106 /* -------------------------------------------------------------------*/ 2107 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2108 MatGetRow_SeqAIJ, 2109 MatRestoreRow_SeqAIJ, 2110 MatMult_SeqAIJ, 2111 MatMultAdd_SeqAIJ, 2112 MatMultTranspose_SeqAIJ, 2113 MatMultTransposeAdd_SeqAIJ, 2114 MatSolve_SeqAIJ, 2115 MatSolveAdd_SeqAIJ, 2116 MatSolveTranspose_SeqAIJ, 2117 MatSolveTransposeAdd_SeqAIJ, 2118 MatLUFactor_SeqAIJ, 2119 0, 2120 MatRelax_SeqAIJ, 2121 MatTranspose_SeqAIJ, 2122 MatGetInfo_SeqAIJ, 2123 MatEqual_SeqAIJ, 2124 MatGetDiagonal_SeqAIJ, 2125 MatDiagonalScale_SeqAIJ, 2126 MatNorm_SeqAIJ, 2127 0, 2128 MatAssemblyEnd_SeqAIJ, 2129 MatCompress_SeqAIJ, 2130 MatSetOption_SeqAIJ, 2131 MatZeroEntries_SeqAIJ, 2132 MatZeroRows_SeqAIJ, 2133 MatLUFactorSymbolic_SeqAIJ, 2134 MatLUFactorNumeric_SeqAIJ, 2135 MatCholeskyFactorSymbolic_SeqAIJ, 2136 MatCholeskyFactorNumeric_SeqAIJ, 2137 MatSetUpPreallocation_SeqAIJ, 2138 MatILUFactorSymbolic_SeqAIJ, 2139 MatICCFactorSymbolic_SeqAIJ, 2140 MatGetArray_SeqAIJ, 2141 MatRestoreArray_SeqAIJ, 2142 MatDuplicate_SeqAIJ, 2143 0, 2144 0, 2145 MatILUFactor_SeqAIJ, 2146 0, 2147 MatAXPY_SeqAIJ, 2148 MatGetSubMatrices_SeqAIJ, 2149 MatIncreaseOverlap_SeqAIJ, 2150 MatGetValues_SeqAIJ, 2151 MatCopy_SeqAIJ, 2152 MatPrintHelp_SeqAIJ, 2153 MatScale_SeqAIJ, 2154 0, 2155 0, 2156 MatILUDTFactor_SeqAIJ, 2157 MatGetBlockSize_SeqAIJ, 2158 MatGetRowIJ_SeqAIJ, 2159 MatRestoreRowIJ_SeqAIJ, 2160 MatGetColumnIJ_SeqAIJ, 2161 MatRestoreColumnIJ_SeqAIJ, 2162 MatFDColoringCreate_SeqAIJ, 2163 0, 2164 0, 2165 MatPermute_SeqAIJ, 2166 0, 2167 0, 2168 MatDestroy_SeqAIJ, 2169 MatView_SeqAIJ, 2170 MatGetPetscMaps_Petsc, 2171 0, 2172 0, 2173 0, 2174 0, 2175 0, 2176 0, 2177 0, 2178 0, 2179 MatSetColoring_SeqAIJ, 2180 MatSetValuesAdic_SeqAIJ, 2181 MatSetValuesAdifor_SeqAIJ, 2182 MatFDColoringApply_SeqAIJ}; 2183 2184 EXTERN_C_BEGIN 2185 #undef __FUNCT__ 2186 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2187 2188 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 2189 { 2190 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2191 int i,nz,n; 2192 2193 PetscFunctionBegin; 2194 2195 nz = aij->maxnz; 2196 n = mat->n; 2197 for (i=0; i<nz; i++) { 2198 aij->j[i] = indices[i]; 2199 } 2200 aij->nz = nz; 2201 for (i=0; i<n; i++) { 2202 aij->ilen[i] = aij->imax[i]; 2203 } 2204 2205 PetscFunctionReturn(0); 2206 } 2207 EXTERN_C_END 2208 2209 #undef __FUNCT__ 2210 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2211 /*@ 2212 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2213 in the matrix. 2214 2215 Input Parameters: 2216 + mat - the SeqAIJ matrix 2217 - indices - the column indices 2218 2219 Level: advanced 2220 2221 Notes: 2222 This can be called if you have precomputed the nonzero structure of the 2223 matrix and want to provide it to the matrix object to improve the performance 2224 of the MatSetValues() operation. 2225 2226 You MUST have set the correct numbers of nonzeros per row in the call to 2227 MatCreateSeqAIJ(). 2228 2229 MUST be called before any calls to MatSetValues(); 2230 2231 The indices should start with zero, not one. 2232 2233 @*/ 2234 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 2235 { 2236 int ierr,(*f)(Mat,int *); 2237 2238 PetscFunctionBegin; 2239 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2240 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2241 if (f) { 2242 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2243 } else { 2244 SETERRQ(1,"Wrong type of matrix to set column indices"); 2245 } 2246 PetscFunctionReturn(0); 2247 } 2248 2249 /* ----------------------------------------------------------------------------------------*/ 2250 2251 EXTERN_C_BEGIN 2252 #undef __FUNCT__ 2253 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2254 int MatStoreValues_SeqAIJ(Mat mat) 2255 { 2256 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2257 size_t nz = aij->i[mat->m],ierr; 2258 2259 PetscFunctionBegin; 2260 if (aij->nonew != 1) { 2261 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2262 } 2263 2264 /* allocate space for values if not already there */ 2265 if (!aij->saved_values) { 2266 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2267 } 2268 2269 /* copy values over */ 2270 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2271 PetscFunctionReturn(0); 2272 } 2273 EXTERN_C_END 2274 2275 #undef __FUNCT__ 2276 #define __FUNCT__ "MatStoreValues" 2277 /*@ 2278 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2279 example, reuse of the linear part of a Jacobian, while recomputing the 2280 nonlinear portion. 2281 2282 Collect on Mat 2283 2284 Input Parameters: 2285 . mat - the matrix (currently on AIJ matrices support this option) 2286 2287 Level: advanced 2288 2289 Common Usage, with SNESSolve(): 2290 $ Create Jacobian matrix 2291 $ Set linear terms into matrix 2292 $ Apply boundary conditions to matrix, at this time matrix must have 2293 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2294 $ boundary conditions again will not change the nonzero structure 2295 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2296 $ ierr = MatStoreValues(mat); 2297 $ Call SNESSetJacobian() with matrix 2298 $ In your Jacobian routine 2299 $ ierr = MatRetrieveValues(mat); 2300 $ Set nonlinear terms in matrix 2301 2302 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2303 $ // build linear portion of Jacobian 2304 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2305 $ ierr = MatStoreValues(mat); 2306 $ loop over nonlinear iterations 2307 $ ierr = MatRetrieveValues(mat); 2308 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2309 $ // call MatAssemblyBegin/End() on matrix 2310 $ Solve linear system with Jacobian 2311 $ endloop 2312 2313 Notes: 2314 Matrix must already be assemblied before calling this routine 2315 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2316 calling this routine. 2317 2318 .seealso: MatRetrieveValues() 2319 2320 @*/ 2321 int MatStoreValues(Mat mat) 2322 { 2323 int ierr,(*f)(Mat); 2324 2325 PetscFunctionBegin; 2326 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2327 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2328 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2329 2330 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2331 if (f) { 2332 ierr = (*f)(mat);CHKERRQ(ierr); 2333 } else { 2334 SETERRQ(1,"Wrong type of matrix to store values"); 2335 } 2336 PetscFunctionReturn(0); 2337 } 2338 2339 EXTERN_C_BEGIN 2340 #undef __FUNCT__ 2341 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2342 int MatRetrieveValues_SeqAIJ(Mat mat) 2343 { 2344 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2345 int nz = aij->i[mat->m],ierr; 2346 2347 PetscFunctionBegin; 2348 if (aij->nonew != 1) { 2349 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2350 } 2351 if (!aij->saved_values) { 2352 SETERRQ(1,"Must call MatStoreValues(A);first"); 2353 } 2354 2355 /* copy values over */ 2356 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2357 PetscFunctionReturn(0); 2358 } 2359 EXTERN_C_END 2360 2361 #undef __FUNCT__ 2362 #define __FUNCT__ "MatRetrieveValues" 2363 /*@ 2364 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2365 example, reuse of the linear part of a Jacobian, while recomputing the 2366 nonlinear portion. 2367 2368 Collect on Mat 2369 2370 Input Parameters: 2371 . mat - the matrix (currently on AIJ matrices support this option) 2372 2373 Level: advanced 2374 2375 .seealso: MatStoreValues() 2376 2377 @*/ 2378 int MatRetrieveValues(Mat mat) 2379 { 2380 int ierr,(*f)(Mat); 2381 2382 PetscFunctionBegin; 2383 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2384 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2385 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2386 2387 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2388 if (f) { 2389 ierr = (*f)(mat);CHKERRQ(ierr); 2390 } else { 2391 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2392 } 2393 PetscFunctionReturn(0); 2394 } 2395 2396 /* 2397 This allows SeqAIJ matrices to be passed to the matlab engine 2398 */ 2399 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2400 #include "engine.h" /* Matlab include file */ 2401 #include "mex.h" /* Matlab include file */ 2402 EXTERN_C_BEGIN 2403 #undef __FUNCT__ 2404 #define __FUNCT__ "MatMatlabEnginePut_SeqAIJ" 2405 int MatMatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine) 2406 { 2407 int ierr,i,*ai,*aj; 2408 Mat B = (Mat)obj; 2409 mxArray *mat; 2410 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2411 2412 PetscFunctionBegin; 2413 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2414 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2415 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2416 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2417 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2418 2419 /* Matlab indices start at 0 for sparse (what a surprise) */ 2420 2421 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2422 mxSetName(mat,obj->name); 2423 engPutArray((Engine *)mengine,mat); 2424 PetscFunctionReturn(0); 2425 } 2426 EXTERN_C_END 2427 2428 EXTERN_C_BEGIN 2429 #undef __FUNCT__ 2430 #define __FUNCT__ "MatMatlabEngineGet_SeqAIJ" 2431 int MatMatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine) 2432 { 2433 int ierr,ii; 2434 Mat mat = (Mat)obj; 2435 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2436 mxArray *mmat; 2437 2438 PetscFunctionBegin; 2439 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2440 2441 mmat = engGetArray((Engine *)mengine,obj->name); 2442 2443 aij->nz = (mxGetJc(mmat))[mat->m]; 2444 ierr = PetscMalloc(((size_t) aij->nz)*(sizeof(int)+sizeof(PetscScalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 2445 aij->j = (int*)(aij->a + aij->nz); 2446 aij->i = aij->j + aij->nz; 2447 aij->singlemalloc = PETSC_TRUE; 2448 aij->freedata = PETSC_TRUE; 2449 2450 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2451 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2452 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 2453 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 2454 2455 for (ii=0; ii<mat->m; ii++) { 2456 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 2457 } 2458 2459 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2460 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2461 2462 PetscFunctionReturn(0); 2463 } 2464 EXTERN_C_END 2465 #endif 2466 2467 /* --------------------------------------------------------------------------------*/ 2468 #undef __FUNCT__ 2469 #define __FUNCT__ "MatCreateSeqAIJ" 2470 /*@C 2471 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2472 (the default parallel PETSc format). For good matrix assembly performance 2473 the user should preallocate the matrix storage by setting the parameter nz 2474 (or the array nnz). By setting these parameters accurately, performance 2475 during matrix assembly can be increased by more than a factor of 50. 2476 2477 Collective on MPI_Comm 2478 2479 Input Parameters: 2480 + comm - MPI communicator, set to PETSC_COMM_SELF 2481 . m - number of rows 2482 . n - number of columns 2483 . nz - number of nonzeros per row (same for all rows) 2484 - nnz - array containing the number of nonzeros in the various rows 2485 (possibly different for each row) or PETSC_NULL 2486 2487 Output Parameter: 2488 . A - the matrix 2489 2490 Notes: 2491 The AIJ format (also called the Yale sparse matrix format or 2492 compressed row storage), is fully compatible with standard Fortran 77 2493 storage. That is, the stored row and column indices can begin at 2494 either one (as in Fortran) or zero. See the users' manual for details. 2495 2496 Specify the preallocated storage with either nz or nnz (not both). 2497 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2498 allocation. For large problems you MUST preallocate memory or you 2499 will get TERRIBLE performance, see the users' manual chapter on matrices. 2500 2501 By default, this format uses inodes (identical nodes) when possible, to 2502 improve numerical efficiency of matrix-vector products and solves. We 2503 search for consecutive rows with the same nonzero structure, thereby 2504 reusing matrix information to achieve increased efficiency. 2505 2506 Options Database Keys: 2507 + -mat_aij_no_inode - Do not use inodes 2508 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2509 - -mat_aij_oneindex - Internally use indexing starting at 1 2510 rather than 0. Note that when calling MatSetValues(), 2511 the user still MUST index entries starting at 0! 2512 2513 Level: intermediate 2514 2515 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2516 2517 @*/ 2518 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz,Mat *A) 2519 { 2520 int ierr; 2521 2522 PetscFunctionBegin; 2523 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2524 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2525 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2526 PetscFunctionReturn(0); 2527 } 2528 2529 #define SKIP_ALLOCATION -4 2530 2531 #undef __FUNCT__ 2532 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2533 /*@C 2534 MatSeqAIJSetPreallocation - For good matrix assembly performance 2535 the user should preallocate the matrix storage by setting the parameter nz 2536 (or the array nnz). By setting these parameters accurately, performance 2537 during matrix assembly can be increased by more than a factor of 50. 2538 2539 Collective on MPI_Comm 2540 2541 Input Parameters: 2542 + comm - MPI communicator, set to PETSC_COMM_SELF 2543 . m - number of rows 2544 . n - number of columns 2545 . nz - number of nonzeros per row (same for all rows) 2546 - nnz - array containing the number of nonzeros in the various rows 2547 (possibly different for each row) or PETSC_NULL 2548 2549 Output Parameter: 2550 . A - the matrix 2551 2552 Notes: 2553 The AIJ format (also called the Yale sparse matrix format or 2554 compressed row storage), is fully compatible with standard Fortran 77 2555 storage. That is, the stored row and column indices can begin at 2556 either one (as in Fortran) or zero. See the users' manual for details. 2557 2558 Specify the preallocated storage with either nz or nnz (not both). 2559 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2560 allocation. For large problems you MUST preallocate memory or you 2561 will get TERRIBLE performance, see the users' manual chapter on matrices. 2562 2563 By default, this format uses inodes (identical nodes) when possible, to 2564 improve numerical efficiency of matrix-vector products and solves. We 2565 search for consecutive rows with the same nonzero structure, thereby 2566 reusing matrix information to achieve increased efficiency. 2567 2568 Options Database Keys: 2569 + -mat_aij_no_inode - Do not use inodes 2570 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2571 - -mat_aij_oneindex - Internally use indexing starting at 1 2572 rather than 0. Note that when calling MatSetValues(), 2573 the user still MUST index entries starting at 0! 2574 2575 Level: intermediate 2576 2577 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2578 2579 @*/ 2580 int MatSeqAIJSetPreallocation(Mat B,int nz,int *nnz) 2581 { 2582 int ierr,(*f)(Mat,int,int*); 2583 2584 PetscFunctionBegin; 2585 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2586 if (f) { 2587 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2588 } 2589 PetscFunctionReturn(0); 2590 } 2591 2592 EXTERN_C_BEGIN 2593 #undef __FUNCT__ 2594 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2595 int MatSeqAIJSetPreallocation_SeqAIJ(Mat B,int nz,int *nnz) 2596 { 2597 Mat_SeqAIJ *b; 2598 size_t len = 0; 2599 PetscTruth flg2,skipallocation = PETSC_FALSE; 2600 int i,ierr; 2601 2602 PetscFunctionBegin; 2603 2604 if (nz == SKIP_ALLOCATION) { 2605 skipallocation = PETSC_TRUE; 2606 nz = 0; 2607 } 2608 2609 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2610 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2611 if (nnz) { 2612 for (i=0; i<B->m; i++) { 2613 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2614 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); 2615 } 2616 } 2617 2618 B->preallocated = PETSC_TRUE; 2619 b = (Mat_SeqAIJ*)B->data; 2620 2621 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2622 if (!nnz) { 2623 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2624 else if (nz <= 0) nz = 1; 2625 for (i=0; i<B->m; i++) b->imax[i] = nz; 2626 nz = nz*B->m; 2627 } else { 2628 nz = 0; 2629 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2630 } 2631 2632 if (!skipallocation) { 2633 /* allocate the matrix space */ 2634 len = ((size_t) nz)*(sizeof(int) + sizeof(PetscScalar)) + (B->m+1)*sizeof(int); 2635 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2636 b->j = (int*)(b->a + nz); 2637 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2638 b->i = b->j + nz; 2639 b->i[0] = 0; 2640 for (i=1; i<B->m+1; i++) { 2641 b->i[i] = b->i[i-1] + b->imax[i-1]; 2642 } 2643 b->singlemalloc = PETSC_TRUE; 2644 b->freedata = PETSC_TRUE; 2645 } else { 2646 b->freedata = PETSC_FALSE; 2647 } 2648 2649 /* b->ilen will count nonzeros in each row so far. */ 2650 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2651 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2652 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2653 2654 b->nz = 0; 2655 b->maxnz = nz; 2656 B->info.nz_unneeded = (double)b->maxnz; 2657 PetscFunctionReturn(0); 2658 } 2659 EXTERN_C_END 2660 2661 EXTERN int RegisterApplyPtAPRoutines_Private(Mat); 2662 2663 EXTERN_C_BEGIN 2664 extern int MatConvert_SeqAIJ_SeqSBAIJ(Mat,MatType,Mat*); 2665 extern int MatConvert_SeqAIJ_SeqBAIJ(Mat,MatType,Mat*); 2666 EXTERN_C_END 2667 2668 EXTERN_C_BEGIN 2669 #undef __FUNCT__ 2670 #define __FUNCT__ "MatCreate_SeqAIJ" 2671 int MatCreate_SeqAIJ(Mat B) 2672 { 2673 Mat_SeqAIJ *b; 2674 int ierr,size; 2675 PetscTruth flg; 2676 2677 PetscFunctionBegin; 2678 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2679 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2680 2681 B->m = B->M = PetscMax(B->m,B->M); 2682 B->n = B->N = PetscMax(B->n,B->N); 2683 2684 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2685 B->data = (void*)b; 2686 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2687 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2688 B->factor = 0; 2689 B->lupivotthreshold = 1.0; 2690 B->mapping = 0; 2691 ierr = PetscOptionsGetReal(B->prefix,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2692 ierr = PetscOptionsHasName(B->prefix,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2693 b->row = 0; 2694 b->col = 0; 2695 b->icol = 0; 2696 b->reallocs = 0; 2697 2698 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2699 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2700 2701 b->sorted = PETSC_FALSE; 2702 b->ignorezeroentries = PETSC_FALSE; 2703 b->roworiented = PETSC_TRUE; 2704 b->nonew = 0; 2705 b->diag = 0; 2706 b->solve_work = 0; 2707 B->spptr = 0; 2708 b->inode.use = PETSC_TRUE; 2709 b->inode.node_count = 0; 2710 b->inode.size = 0; 2711 b->inode.limit = 5; 2712 b->inode.max_limit = 5; 2713 b->saved_values = 0; 2714 b->idiag = 0; 2715 b->ssor = 0; 2716 b->keepzeroedrows = PETSC_FALSE; 2717 b->xtoy = 0; 2718 b->XtoY = 0; 2719 2720 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2721 2722 #if defined(PETSC_HAVE_SUPERLU) 2723 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_superlu",&flg);CHKERRQ(ierr); 2724 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B);CHKERRQ(ierr); } 2725 #endif 2726 2727 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_essl",&flg);CHKERRQ(ierr); 2728 if (flg) { ierr = MatUseEssl_SeqAIJ(B);CHKERRQ(ierr); } 2729 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_lusol",&flg);CHKERRQ(ierr); 2730 if (flg) { ierr = MatUseLUSOL_SeqAIJ(B);CHKERRQ(ierr); } 2731 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_matlab",&flg);CHKERRQ(ierr); 2732 if (flg) {ierr = MatUseMatlab_SeqAIJ(B);CHKERRQ(ierr);} 2733 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_dxml",&flg);CHKERRQ(ierr); 2734 if (flg){ 2735 ierr = MatUseDXML_SeqAIJ(B);CHKERRQ(ierr); 2736 } 2737 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2738 "MatSeqAIJSetColumnIndices_SeqAIJ", 2739 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2740 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2741 "MatStoreValues_SeqAIJ", 2742 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2743 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2744 "MatRetrieveValues_SeqAIJ", 2745 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2746 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2747 "MatConvert_SeqAIJ_SeqSBAIJ", 2748 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2749 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2750 "MatConvert_SeqAIJ_SeqBAIJ", 2751 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2752 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsSymmetric_C", 2753 "MatIsSymmetric_SeqAIJ", 2754 MatIsSymmetric_SeqAIJ);CHKERRQ(ierr); 2755 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2756 "MatSeqAIJSetPreallocation_SeqAIJ", 2757 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2758 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2759 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2760 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2761 #endif 2762 ierr = RegisterApplyPtAPRoutines_Private(B);CHKERRQ(ierr); 2763 PetscFunctionReturn(0); 2764 } 2765 EXTERN_C_END 2766 2767 #undef __FUNCT__ 2768 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2769 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2770 { 2771 Mat C; 2772 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2773 int i,m = A->m,ierr; 2774 size_t len; 2775 2776 PetscFunctionBegin; 2777 *B = 0; 2778 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2779 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 2780 c = (Mat_SeqAIJ*)C->data; 2781 2782 C->factor = A->factor; 2783 c->row = 0; 2784 c->col = 0; 2785 c->icol = 0; 2786 c->keepzeroedrows = a->keepzeroedrows; 2787 C->assembled = PETSC_TRUE; 2788 2789 C->M = A->m; 2790 C->N = A->n; 2791 2792 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2793 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2794 for (i=0; i<m; i++) { 2795 c->imax[i] = a->imax[i]; 2796 c->ilen[i] = a->ilen[i]; 2797 } 2798 2799 /* allocate the matrix space */ 2800 c->singlemalloc = PETSC_TRUE; 2801 len = ((size_t) (m+1))*sizeof(int)+(a->i[m])*(sizeof(PetscScalar)+sizeof(int)); 2802 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2803 c->j = (int*)(c->a + a->i[m] ); 2804 c->i = c->j + a->i[m]; 2805 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2806 if (m > 0) { 2807 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(int));CHKERRQ(ierr); 2808 if (cpvalues == MAT_COPY_VALUES) { 2809 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2810 } else { 2811 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2812 } 2813 } 2814 2815 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2816 c->sorted = a->sorted; 2817 c->roworiented = a->roworiented; 2818 c->nonew = a->nonew; 2819 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2820 c->saved_values = 0; 2821 c->idiag = 0; 2822 c->ssor = 0; 2823 c->ignorezeroentries = a->ignorezeroentries; 2824 c->freedata = PETSC_TRUE; 2825 2826 if (a->diag) { 2827 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2828 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2829 for (i=0; i<m; i++) { 2830 c->diag[i] = a->diag[i]; 2831 } 2832 } else c->diag = 0; 2833 c->inode.use = a->inode.use; 2834 c->inode.limit = a->inode.limit; 2835 c->inode.max_limit = a->inode.max_limit; 2836 if (a->inode.size){ 2837 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2838 c->inode.node_count = a->inode.node_count; 2839 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2840 } else { 2841 c->inode.size = 0; 2842 c->inode.node_count = 0; 2843 } 2844 c->nz = a->nz; 2845 c->maxnz = a->maxnz; 2846 c->solve_work = 0; 2847 C->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2848 C->preallocated = PETSC_TRUE; 2849 2850 *B = C; 2851 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2852 PetscFunctionReturn(0); 2853 } 2854 2855 EXTERN_C_BEGIN 2856 #undef __FUNCT__ 2857 #define __FUNCT__ "MatLoad_SeqAIJ" 2858 int MatLoad_SeqAIJ(PetscViewer viewer,MatType type,Mat *A) 2859 { 2860 Mat_SeqAIJ *a; 2861 Mat B; 2862 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N; 2863 MPI_Comm comm; 2864 #if defined(PETSC_HAVE_SPOOLES) || defined(PETSC_HAVE_SUPERLU) || defined(PETSC_HAVE_SUPERLUDIST) || defined(PETSC_HAVE_UMFPACK) || defined(PETSC_HAVE_MUMPS) 2865 PetscTruth flag; 2866 #endif 2867 2868 PetscFunctionBegin; 2869 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2870 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2871 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2872 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2873 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2874 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2875 M = header[1]; N = header[2]; nz = header[3]; 2876 2877 if (nz < 0) { 2878 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2879 } 2880 2881 /* read in row lengths */ 2882 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2883 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2884 2885 /* create our matrix */ 2886 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A);CHKERRQ(ierr); 2887 B = *A; 2888 a = (Mat_SeqAIJ*)B->data; 2889 2890 /* read in column indices and adjust for Fortran indexing*/ 2891 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2892 2893 /* read in nonzero values */ 2894 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2895 2896 /* set matrix "i" values */ 2897 a->i[0] = 0; 2898 for (i=1; i<= M; i++) { 2899 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2900 a->ilen[i-1] = rowlengths[i-1]; 2901 } 2902 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2903 2904 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2905 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2906 #if defined(PETSC_HAVE_SPOOLES) 2907 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_spooles",&flag);CHKERRQ(ierr); 2908 if (flag) { ierr = MatUseSpooles_SeqAIJ(B);CHKERRQ(ierr); } 2909 #endif 2910 #if defined(PETSC_HAVE_SUPERLU) 2911 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_superlu",&flag);CHKERRQ(ierr); 2912 if (flag) { ierr = MatUseSuperLU_SeqAIJ(B);CHKERRQ(ierr); } 2913 #endif 2914 #if defined(PETSC_HAVE_SUPERLUDIST) 2915 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_superlu_dist",&flag);CHKERRQ(ierr); 2916 if (flag) { ierr = MatUseSuperLU_DIST_MPIAIJ(B);CHKERRQ(ierr); } 2917 #endif 2918 #if defined(PETSC_HAVE_UMFPACK) 2919 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_umfpack",&flag);CHKERRQ(ierr); 2920 if (flag) { ierr = MatUseUMFPACK_SeqAIJ(B);CHKERRQ(ierr); } 2921 #endif 2922 #if defined(PETSC_HAVE_MUMPS) 2923 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_mumps",&flag);CHKERRQ(ierr); 2924 if (flag) { ierr = MatUseMUMPS_MPIAIJ(B);CHKERRQ(ierr); } 2925 #endif 2926 PetscFunctionReturn(0); 2927 } 2928 EXTERN_C_END 2929 2930 #undef __FUNCT__ 2931 #define __FUNCT__ "MatEqual_SeqAIJ" 2932 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2933 { 2934 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2935 int ierr; 2936 2937 PetscFunctionBegin; 2938 /* If the matrix dimensions are not equal,or no of nonzeros */ 2939 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)) { 2940 *flg = PETSC_FALSE; 2941 PetscFunctionReturn(0); 2942 } 2943 2944 /* if the a->i are the same */ 2945 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2946 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2947 2948 /* if a->j are the same */ 2949 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2950 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2951 2952 /* if a->a are the same */ 2953 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2954 2955 PetscFunctionReturn(0); 2956 2957 } 2958 2959 #undef __FUNCT__ 2960 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2961 /*@C 2962 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2963 provided by the user. 2964 2965 Coolective on MPI_Comm 2966 2967 Input Parameters: 2968 + comm - must be an MPI communicator of size 1 2969 . m - number of rows 2970 . n - number of columns 2971 . i - row indices 2972 . j - column indices 2973 - a - matrix values 2974 2975 Output Parameter: 2976 . mat - the matrix 2977 2978 Level: intermediate 2979 2980 Notes: 2981 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2982 once the matrix is destroyed 2983 2984 You cannot set new nonzero locations into this matrix, that will generate an error. 2985 2986 The i and j indices are 0 based 2987 2988 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2989 2990 @*/ 2991 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,PetscScalar *a,Mat *mat) 2992 { 2993 int ierr,ii; 2994 Mat_SeqAIJ *aij; 2995 2996 PetscFunctionBegin; 2997 ierr = MatCreateSeqAIJ(comm,m,n,SKIP_ALLOCATION,0,mat);CHKERRQ(ierr); 2998 aij = (Mat_SeqAIJ*)(*mat)->data; 2999 3000 if (i[0] != 0) { 3001 SETERRQ(1,"i (row indices) must start with 0"); 3002 } 3003 aij->i = i; 3004 aij->j = j; 3005 aij->a = a; 3006 aij->singlemalloc = PETSC_FALSE; 3007 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3008 aij->freedata = PETSC_FALSE; 3009 3010 for (ii=0; ii<m; ii++) { 3011 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3012 #if defined(PETSC_USE_BOPT_g) 3013 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]); 3014 #endif 3015 } 3016 #if defined(PETSC_USE_BOPT_g) 3017 for (ii=0; ii<aij->i[m]; ii++) { 3018 if (j[ii] < 0) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 3019 if (j[ii] > n - 1) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 3020 } 3021 #endif 3022 3023 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3024 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3025 PetscFunctionReturn(0); 3026 } 3027 3028 #undef __FUNCT__ 3029 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3030 int MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3031 { 3032 int ierr; 3033 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3034 3035 PetscFunctionBegin; 3036 if (coloring->ctype == IS_COLORING_LOCAL) { 3037 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3038 a->coloring = coloring; 3039 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3040 int i,*larray; 3041 ISColoring ocoloring; 3042 ISColoringValue *colors; 3043 3044 /* set coloring for diagonal portion */ 3045 ierr = PetscMalloc((A->n+1)*sizeof(int),&larray);CHKERRQ(ierr); 3046 for (i=0; i<A->n; i++) { 3047 larray[i] = i; 3048 } 3049 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3050 ierr = PetscMalloc((A->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3051 for (i=0; i<A->n; i++) { 3052 colors[i] = coloring->colors[larray[i]]; 3053 } 3054 ierr = PetscFree(larray);CHKERRQ(ierr); 3055 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 3056 a->coloring = ocoloring; 3057 } 3058 PetscFunctionReturn(0); 3059 } 3060 3061 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 3062 EXTERN_C_BEGIN 3063 #include "adic/ad_utils.h" 3064 EXTERN_C_END 3065 3066 #undef __FUNCT__ 3067 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3068 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3069 { 3070 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3071 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3072 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3073 ISColoringValue *color; 3074 3075 PetscFunctionBegin; 3076 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3077 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3078 color = a->coloring->colors; 3079 /* loop over rows */ 3080 for (i=0; i<m; i++) { 3081 nz = ii[i+1] - ii[i]; 3082 /* loop over columns putting computed value into matrix */ 3083 for (j=0; j<nz; j++) { 3084 *v++ = values[color[*jj++]]; 3085 } 3086 values += nlen; /* jump to next row of derivatives */ 3087 } 3088 PetscFunctionReturn(0); 3089 } 3090 3091 #else 3092 3093 #undef __FUNCT__ 3094 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3095 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3096 { 3097 PetscFunctionBegin; 3098 SETERRQ(1,"PETSc installed without ADIC"); 3099 } 3100 3101 #endif 3102 3103 #undef __FUNCT__ 3104 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3105 int MatSetValuesAdifor_SeqAIJ(Mat A,int nl,void *advalues) 3106 { 3107 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3108 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j; 3109 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3110 ISColoringValue *color; 3111 3112 PetscFunctionBegin; 3113 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3114 color = a->coloring->colors; 3115 /* loop over rows */ 3116 for (i=0; i<m; i++) { 3117 nz = ii[i+1] - ii[i]; 3118 /* loop over columns putting computed value into matrix */ 3119 for (j=0; j<nz; j++) { 3120 *v++ = values[color[*jj++]]; 3121 } 3122 values += nl; /* jump to next row of derivatives */ 3123 } 3124 PetscFunctionReturn(0); 3125 } 3126 3127