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