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