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