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