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