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 PetscScalar *array; 2387 mxArray *mat; 2388 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2389 2390 PetscFunctionBegin; 2391 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2392 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2393 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2394 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2395 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2396 2397 /* Matlab indices start at 0 for sparse (what a surprise) */ 2398 if (aij->indexshift) { 2399 for (i=0; i<B->m+1; i++) { 2400 ai[i]--; 2401 } 2402 for (i=0; i<aij->nz; i++) { 2403 aj[i]--; 2404 } 2405 } 2406 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2407 mxSetName(mat,obj->name); 2408 engPutArray((Engine *)engine,mat); 2409 PetscFunctionReturn(0); 2410 } 2411 EXTERN_C_END 2412 2413 EXTERN_C_BEGIN 2414 #undef __FUNCT__ 2415 #define __FUNCT__ "MatMatlabEngineGet_SeqAIJ" 2416 int MatMatlabEngineGet_SeqAIJ(PetscObject obj,void *engine) 2417 { 2418 int ierr,ii; 2419 Mat mat = (Mat)obj; 2420 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2421 mxArray *mmat; 2422 2423 PetscFunctionBegin; 2424 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2425 aij->indexshift = 0; 2426 2427 mmat = engGetArray((Engine *)engine,obj->name); 2428 2429 aij->nz = (mxGetJc(mmat))[mat->m]; 2430 ierr = PetscMalloc(aij->nz*(sizeof(int)+sizeof(PetscScalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 2431 aij->j = (int*)(aij->a + aij->nz); 2432 aij->i = aij->j + aij->nz; 2433 aij->singlemalloc = PETSC_TRUE; 2434 aij->freedata = PETSC_TRUE; 2435 2436 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2437 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2438 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 2439 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 2440 2441 for (ii=0; ii<mat->m; ii++) { 2442 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 2443 } 2444 2445 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2446 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2447 2448 PetscFunctionReturn(0); 2449 } 2450 EXTERN_C_END 2451 #endif 2452 2453 /* --------------------------------------------------------------------------------*/ 2454 #undef __FUNCT__ 2455 #define __FUNCT__ "MatCreateSeqAIJ" 2456 /*@C 2457 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2458 (the default parallel PETSc format). For good matrix assembly performance 2459 the user should preallocate the matrix storage by setting the parameter nz 2460 (or the array nnz). By setting these parameters accurately, performance 2461 during matrix assembly can be increased by more than a factor of 50. 2462 2463 Collective on MPI_Comm 2464 2465 Input Parameters: 2466 + comm - MPI communicator, set to PETSC_COMM_SELF 2467 . m - number of rows 2468 . n - number of columns 2469 . nz - number of nonzeros per row (same for all rows) 2470 - nnz - array containing the number of nonzeros in the various rows 2471 (possibly different for each row) or PETSC_NULL 2472 2473 Output Parameter: 2474 . A - the matrix 2475 2476 Notes: 2477 The AIJ format (also called the Yale sparse matrix format or 2478 compressed row storage), is fully compatible with standard Fortran 77 2479 storage. That is, the stored row and column indices can begin at 2480 either one (as in Fortran) or zero. See the users' manual for details. 2481 2482 Specify the preallocated storage with either nz or nnz (not both). 2483 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2484 allocation. For large problems you MUST preallocate memory or you 2485 will get TERRIBLE performance, see the users' manual chapter on matrices. 2486 2487 By default, this format uses inodes (identical nodes) when possible, to 2488 improve numerical efficiency of matrix-vector products and solves. We 2489 search for consecutive rows with the same nonzero structure, thereby 2490 reusing matrix information to achieve increased efficiency. 2491 2492 Options Database Keys: 2493 + -mat_aij_no_inode - Do not use inodes 2494 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2495 - -mat_aij_oneindex - Internally use indexing starting at 1 2496 rather than 0. Note that when calling MatSetValues(), 2497 the user still MUST index entries starting at 0! 2498 2499 Level: intermediate 2500 2501 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2502 2503 @*/ 2504 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz,Mat *A) 2505 { 2506 int ierr; 2507 2508 PetscFunctionBegin; 2509 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2510 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2511 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2512 PetscFunctionReturn(0); 2513 } 2514 2515 #define SKIP_ALLOCATION -4 2516 2517 #undef __FUNCT__ 2518 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2519 /*@C 2520 MatSeqAIJSetPreallocation - For good matrix assembly performance 2521 the user should preallocate the matrix storage by setting the parameter nz 2522 (or the array nnz). By setting these parameters accurately, performance 2523 during matrix assembly can be increased by more than a factor of 50. 2524 2525 Collective on MPI_Comm 2526 2527 Input Parameters: 2528 + comm - MPI communicator, set to PETSC_COMM_SELF 2529 . m - number of rows 2530 . n - number of columns 2531 . nz - number of nonzeros per row (same for all rows) 2532 - nnz - array containing the number of nonzeros in the various rows 2533 (possibly different for each row) or PETSC_NULL 2534 2535 Output Parameter: 2536 . A - the matrix 2537 2538 Notes: 2539 The AIJ format (also called the Yale sparse matrix format or 2540 compressed row storage), is fully compatible with standard Fortran 77 2541 storage. That is, the stored row and column indices can begin at 2542 either one (as in Fortran) or zero. See the users' manual for details. 2543 2544 Specify the preallocated storage with either nz or nnz (not both). 2545 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2546 allocation. For large problems you MUST preallocate memory or you 2547 will get TERRIBLE performance, see the users' manual chapter on matrices. 2548 2549 By default, this format uses inodes (identical nodes) when possible, to 2550 improve numerical efficiency of matrix-vector products and solves. We 2551 search for consecutive rows with the same nonzero structure, thereby 2552 reusing matrix information to achieve increased efficiency. 2553 2554 Options Database Keys: 2555 + -mat_aij_no_inode - Do not use inodes 2556 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2557 - -mat_aij_oneindex - Internally use indexing starting at 1 2558 rather than 0. Note that when calling MatSetValues(), 2559 the user still MUST index entries starting at 0! 2560 2561 Level: intermediate 2562 2563 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2564 2565 @*/ 2566 int MatSeqAIJSetPreallocation(Mat B,int nz,int *nnz) 2567 { 2568 Mat_SeqAIJ *b; 2569 int i,len=0,ierr; 2570 PetscTruth flg2,skipallocation = PETSC_FALSE; 2571 2572 PetscFunctionBegin; 2573 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flg2);CHKERRQ(ierr); 2574 if (!flg2) PetscFunctionReturn(0); 2575 2576 if (nz == SKIP_ALLOCATION) { 2577 skipallocation = PETSC_TRUE; 2578 nz = 0; 2579 } 2580 2581 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2582 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2583 if (nnz) { 2584 for (i=0; i<B->m; i++) { 2585 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2586 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); 2587 } 2588 } 2589 2590 B->preallocated = PETSC_TRUE; 2591 b = (Mat_SeqAIJ*)B->data; 2592 2593 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2594 if (!nnz) { 2595 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2596 else if (nz <= 0) nz = 1; 2597 for (i=0; i<B->m; i++) b->imax[i] = nz; 2598 nz = nz*B->m; 2599 } else { 2600 nz = 0; 2601 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2602 } 2603 2604 if (!skipallocation) { 2605 /* allocate the matrix space */ 2606 len = nz*(sizeof(int) + sizeof(PetscScalar)) + (B->m+1)*sizeof(int); 2607 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2608 b->j = (int*)(b->a + nz); 2609 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2610 b->i = b->j + nz; 2611 b->i[0] = -b->indexshift; 2612 for (i=1; i<B->m+1; i++) { 2613 b->i[i] = b->i[i-1] + b->imax[i-1]; 2614 } 2615 b->singlemalloc = PETSC_TRUE; 2616 b->freedata = PETSC_TRUE; 2617 } else { 2618 b->freedata = PETSC_FALSE; 2619 } 2620 2621 /* b->ilen will count nonzeros in each row so far. */ 2622 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2623 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2624 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2625 2626 b->nz = 0; 2627 b->maxnz = nz; 2628 B->info.nz_unneeded = (double)b->maxnz; 2629 PetscFunctionReturn(0); 2630 } 2631 2632 EXTERN int RegisterApplyPtAPRoutines_Private(Mat); 2633 2634 EXTERN_C_BEGIN 2635 #undef __FUNCT__ 2636 #define __FUNCT__ "MatCreate_SeqAIJ" 2637 int MatCreate_SeqAIJ(Mat B) 2638 { 2639 Mat_SeqAIJ *b; 2640 int ierr,size; 2641 PetscTruth flg; 2642 2643 PetscFunctionBegin; 2644 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2645 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2646 2647 B->m = B->M = PetscMax(B->m,B->M); 2648 B->n = B->N = PetscMax(B->n,B->N); 2649 2650 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2651 B->data = (void*)b; 2652 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2653 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2654 B->factor = 0; 2655 B->lupivotthreshold = 1.0; 2656 B->mapping = 0; 2657 ierr = PetscOptionsGetReal(B->prefix,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2658 ierr = PetscOptionsHasName(B->prefix,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2659 b->row = 0; 2660 b->col = 0; 2661 b->icol = 0; 2662 b->indexshift = 0; 2663 b->reallocs = 0; 2664 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_oneindex",&flg);CHKERRQ(ierr); 2665 if (flg) b->indexshift = -1; 2666 2667 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2668 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2669 2670 b->sorted = PETSC_FALSE; 2671 b->ignorezeroentries = PETSC_FALSE; 2672 b->roworiented = PETSC_TRUE; 2673 b->nonew = 0; 2674 b->diag = 0; 2675 b->solve_work = 0; 2676 B->spptr = 0; 2677 b->inode.use = PETSC_TRUE; 2678 b->inode.node_count = 0; 2679 b->inode.size = 0; 2680 b->inode.limit = 5; 2681 b->inode.max_limit = 5; 2682 b->saved_values = 0; 2683 b->idiag = 0; 2684 b->ssor = 0; 2685 b->keepzeroedrows = PETSC_FALSE; 2686 2687 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2688 2689 #if defined(PETSC_HAVE_SUPERLU) 2690 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_superlu",&flg);CHKERRQ(ierr); 2691 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B);CHKERRQ(ierr); } 2692 #endif 2693 2694 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_essl",&flg);CHKERRQ(ierr); 2695 if (flg) { ierr = MatUseEssl_SeqAIJ(B);CHKERRQ(ierr); } 2696 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_lusol",&flg);CHKERRQ(ierr); 2697 if (flg) { ierr = MatUseLUSOL_SeqAIJ(B);CHKERRQ(ierr); } 2698 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_matlab",&flg);CHKERRQ(ierr); 2699 if (flg) {ierr = MatUseMatlab_SeqAIJ(B);CHKERRQ(ierr);} 2700 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_dxml",&flg);CHKERRQ(ierr); 2701 if (flg) { 2702 if (!b->indexshift) SETERRQ(PETSC_ERR_LIB,"need -mat_aij_oneindex with -mat_aij_dxml"); 2703 ierr = MatUseDXML_SeqAIJ(B);CHKERRQ(ierr); 2704 } 2705 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2706 "MatSeqAIJSetColumnIndices_SeqAIJ", 2707 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2708 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2709 "MatStoreValues_SeqAIJ", 2710 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2711 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2712 "MatRetrieveValues_SeqAIJ", 2713 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2714 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2715 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2716 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2717 #endif 2718 ierr = RegisterApplyPtAPRoutines_Private(B);CHKERRQ(ierr); 2719 PetscFunctionReturn(0); 2720 } 2721 EXTERN_C_END 2722 2723 #undef __FUNCT__ 2724 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2725 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2726 { 2727 Mat C; 2728 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2729 int i,len,m = A->m,shift = a->indexshift,ierr; 2730 2731 PetscFunctionBegin; 2732 *B = 0; 2733 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2734 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 2735 c = (Mat_SeqAIJ*)C->data; 2736 2737 C->factor = A->factor; 2738 c->row = 0; 2739 c->col = 0; 2740 c->icol = 0; 2741 c->indexshift = shift; 2742 c->keepzeroedrows = a->keepzeroedrows; 2743 C->assembled = PETSC_TRUE; 2744 2745 C->M = A->m; 2746 C->N = A->n; 2747 2748 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2749 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2750 for (i=0; i<m; i++) { 2751 c->imax[i] = a->imax[i]; 2752 c->ilen[i] = a->ilen[i]; 2753 } 2754 2755 /* allocate the matrix space */ 2756 c->singlemalloc = PETSC_TRUE; 2757 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(PetscScalar)+sizeof(int)); 2758 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2759 c->j = (int*)(c->a + a->i[m] + shift); 2760 c->i = c->j + a->i[m] + shift; 2761 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2762 if (m > 0) { 2763 ierr = PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int));CHKERRQ(ierr); 2764 if (cpvalues == MAT_COPY_VALUES) { 2765 ierr = PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(PetscScalar));CHKERRQ(ierr); 2766 } else { 2767 ierr = PetscMemzero(c->a,(a->i[m]+shift)*sizeof(PetscScalar));CHKERRQ(ierr); 2768 } 2769 } 2770 2771 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2772 c->sorted = a->sorted; 2773 c->roworiented = a->roworiented; 2774 c->nonew = a->nonew; 2775 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2776 c->saved_values = 0; 2777 c->idiag = 0; 2778 c->ssor = 0; 2779 c->ignorezeroentries = a->ignorezeroentries; 2780 c->freedata = PETSC_TRUE; 2781 2782 if (a->diag) { 2783 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2784 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2785 for (i=0; i<m; i++) { 2786 c->diag[i] = a->diag[i]; 2787 } 2788 } else c->diag = 0; 2789 c->inode.use = a->inode.use; 2790 c->inode.limit = a->inode.limit; 2791 c->inode.max_limit = a->inode.max_limit; 2792 if (a->inode.size){ 2793 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2794 c->inode.node_count = a->inode.node_count; 2795 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2796 } else { 2797 c->inode.size = 0; 2798 c->inode.node_count = 0; 2799 } 2800 c->nz = a->nz; 2801 c->maxnz = a->maxnz; 2802 c->solve_work = 0; 2803 C->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2804 C->preallocated = PETSC_TRUE; 2805 2806 *B = C; 2807 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2808 PetscFunctionReturn(0); 2809 } 2810 2811 EXTERN_C_BEGIN 2812 #undef __FUNCT__ 2813 #define __FUNCT__ "MatLoad_SeqAIJ" 2814 int MatLoad_SeqAIJ(PetscViewer viewer,MatType type,Mat *A) 2815 { 2816 Mat_SeqAIJ *a; 2817 Mat B; 2818 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N,shift; 2819 MPI_Comm comm; 2820 2821 PetscFunctionBegin; 2822 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2823 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2824 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2825 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2826 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2827 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2828 M = header[1]; N = header[2]; nz = header[3]; 2829 2830 if (nz < 0) { 2831 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2832 } 2833 2834 /* read in row lengths */ 2835 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2836 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2837 2838 /* create our matrix */ 2839 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A);CHKERRQ(ierr); 2840 B = *A; 2841 a = (Mat_SeqAIJ*)B->data; 2842 shift = a->indexshift; 2843 2844 /* read in column indices and adjust for Fortran indexing*/ 2845 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2846 if (shift) { 2847 for (i=0; i<nz; i++) { 2848 a->j[i] += 1; 2849 } 2850 } 2851 2852 /* read in nonzero values */ 2853 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2854 2855 /* set matrix "i" values */ 2856 a->i[0] = -shift; 2857 for (i=1; i<= M; i++) { 2858 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2859 a->ilen[i-1] = rowlengths[i-1]; 2860 } 2861 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2862 2863 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2864 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2865 PetscFunctionReturn(0); 2866 } 2867 EXTERN_C_END 2868 2869 #undef __FUNCT__ 2870 #define __FUNCT__ "MatEqual_SeqAIJ" 2871 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2872 { 2873 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2874 int ierr; 2875 PetscTruth flag; 2876 2877 PetscFunctionBegin; 2878 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flag);CHKERRQ(ierr); 2879 if (!flag) SETERRQ(PETSC_ERR_ARG_INCOMP,"Matrices must be same type"); 2880 2881 /* If the matrix dimensions are not equal,or no of nonzeros or shift */ 2882 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)|| (a->indexshift != b->indexshift)) { 2883 *flg = PETSC_FALSE; 2884 PetscFunctionReturn(0); 2885 } 2886 2887 /* if the a->i are the same */ 2888 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2889 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2890 2891 /* if a->j are the same */ 2892 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2893 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2894 2895 /* if a->a are the same */ 2896 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2897 2898 PetscFunctionReturn(0); 2899 2900 } 2901 2902 #undef __FUNCT__ 2903 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2904 /*@C 2905 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2906 provided by the user. 2907 2908 Coolective on MPI_Comm 2909 2910 Input Parameters: 2911 + comm - must be an MPI communicator of size 1 2912 . m - number of rows 2913 . n - number of columns 2914 . i - row indices 2915 . j - column indices 2916 - a - matrix values 2917 2918 Output Parameter: 2919 . mat - the matrix 2920 2921 Level: intermediate 2922 2923 Notes: 2924 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2925 once the matrix is destroyed 2926 2927 You cannot set new nonzero locations into this matrix, that will generate an error. 2928 2929 The i and j indices can be either 0- or 1 based 2930 2931 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2932 2933 @*/ 2934 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,PetscScalar *a,Mat *mat) 2935 { 2936 int ierr,ii; 2937 Mat_SeqAIJ *aij; 2938 2939 PetscFunctionBegin; 2940 ierr = MatCreateSeqAIJ(comm,m,n,SKIP_ALLOCATION,0,mat);CHKERRQ(ierr); 2941 aij = (Mat_SeqAIJ*)(*mat)->data; 2942 2943 if (i[0] == 1) { 2944 aij->indexshift = -1; 2945 } else if (i[0]) { 2946 SETERRQ(1,"i (row indices) do not start with 0 or 1"); 2947 } 2948 aij->i = i; 2949 aij->j = j; 2950 aij->a = a; 2951 aij->singlemalloc = PETSC_FALSE; 2952 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2953 aij->freedata = PETSC_FALSE; 2954 2955 for (ii=0; ii<m; ii++) { 2956 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2957 #if defined(PETSC_USE_BOPT_g) 2958 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]); 2959 #endif 2960 } 2961 #if defined(PETSC_USE_BOPT_g) 2962 for (ii=0; ii<aij->i[m]; ii++) { 2963 if (j[ii] < -aij->indexshift) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 2964 if (j[ii] > n - 1 -aij->indexshift) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 2965 } 2966 #endif 2967 2968 /* changes indices to start at 0 */ 2969 if (i[0]) { 2970 aij->indexshift = 0; 2971 for (ii=0; ii<m; ii++) { 2972 i[ii]--; 2973 } 2974 for (ii=0; ii<i[m]; ii++) { 2975 j[ii]--; 2976 } 2977 } 2978 2979 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2980 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2981 PetscFunctionReturn(0); 2982 } 2983 2984 #undef __FUNCT__ 2985 #define __FUNCT__ "MatSetColoring_SeqAIJ" 2986 int MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 2987 { 2988 int ierr; 2989 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2990 2991 PetscFunctionBegin; 2992 if (coloring->ctype == IS_COLORING_LOCAL) { 2993 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 2994 a->coloring = coloring; 2995 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 2996 int *colors,i,*larray; 2997 ISColoring ocoloring; 2998 2999 /* set coloring for diagonal portion */ 3000 ierr = PetscMalloc((A->n+1)*sizeof(int),&larray);CHKERRQ(ierr); 3001 for (i=0; i<A->n; i++) { 3002 larray[i] = i; 3003 } 3004 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3005 ierr = PetscMalloc((A->n+1)*sizeof(int),&colors);CHKERRQ(ierr); 3006 for (i=0; i<A->n; i++) { 3007 colors[i] = coloring->colors[larray[i]]; 3008 } 3009 ierr = PetscFree(larray);CHKERRQ(ierr); 3010 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 3011 a->coloring = ocoloring; 3012 } 3013 PetscFunctionReturn(0); 3014 } 3015 3016 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 3017 EXTERN_C_BEGIN 3018 #include "adic/ad_utils.h" 3019 EXTERN_C_END 3020 3021 #undef __FUNCT__ 3022 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3023 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3024 { 3025 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3026 int m = A->m,*ii = a->i,*jj = a->j,nz,i,*color,j,nlen; 3027 PetscScalar *v = a->a,*values; 3028 char *cadvalues = (char *)advalues; 3029 3030 PetscFunctionBegin; 3031 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3032 nlen = PetscADGetDerivTypeSize(); 3033 color = a->coloring->colors; 3034 /* loop over rows */ 3035 for (i=0; i<m; i++) { 3036 nz = ii[i+1] - ii[i]; 3037 /* loop over columns putting computed value into matrix */ 3038 values = PetscADGetGradArray(cadvalues); 3039 for (j=0; j<nz; j++) { 3040 *v++ = values[color[*jj++]]; 3041 } 3042 cadvalues += nlen; /* jump to next row of derivatives */ 3043 } 3044 PetscFunctionReturn(0); 3045 } 3046 3047 #else 3048 3049 #undef __FUNCT__ 3050 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3051 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3052 { 3053 PetscFunctionBegin; 3054 SETERRQ(1,"PETSc installed without ADIC"); 3055 } 3056 3057 #endif 3058 3059 #undef __FUNCT__ 3060 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3061 int MatSetValuesAdifor_SeqAIJ(Mat A,int nl,void *advalues) 3062 { 3063 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3064 int m = A->m,*ii = a->i,*jj = a->j,nz,i,*color,j; 3065 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3066 3067 PetscFunctionBegin; 3068 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3069 color = a->coloring->colors; 3070 /* loop over rows */ 3071 for (i=0; i<m; i++) { 3072 nz = ii[i+1] - ii[i]; 3073 /* loop over columns putting computed value into matrix */ 3074 for (j=0; j<nz; j++) { 3075 *v++ = values[color[*jj++]]; 3076 } 3077 values += nl; /* jump to next row of derivatives */ 3078 } 3079 PetscFunctionReturn(0); 3080 } 3081 3082