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