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