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