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