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