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