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