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 0, 2055 MatILUFactorSymbolic_SeqAIJ, 2056 0, 2057 MatGetArray_SeqAIJ, 2058 MatRestoreArray_SeqAIJ, 2059 MatDuplicate_SeqAIJ, 2060 0, 2061 0, 2062 MatILUFactor_SeqAIJ, 2063 0, 2064 0, 2065 MatGetSubMatrices_SeqAIJ, 2066 MatIncreaseOverlap_SeqAIJ, 2067 MatGetValues_SeqAIJ, 2068 MatCopy_SeqAIJ, 2069 MatPrintHelp_SeqAIJ, 2070 MatScale_SeqAIJ, 2071 0, 2072 0, 2073 MatILUDTFactor_SeqAIJ, 2074 MatGetBlockSize_SeqAIJ, 2075 MatGetRowIJ_SeqAIJ, 2076 MatRestoreRowIJ_SeqAIJ, 2077 MatGetColumnIJ_SeqAIJ, 2078 MatRestoreColumnIJ_SeqAIJ, 2079 MatFDColoringCreate_SeqAIJ, 2080 0, 2081 0, 2082 MatPermute_SeqAIJ, 2083 0, 2084 0, 2085 MatDestroy_SeqAIJ, 2086 MatView_SeqAIJ, 2087 MatGetPetscMaps_Petsc, 2088 0, 2089 0, 2090 0, 2091 0, 2092 0, 2093 0, 2094 0, 2095 0, 2096 MatSetColoring_SeqAIJ, 2097 MatSetValuesAdic_SeqAIJ, 2098 MatSetValuesAdifor_SeqAIJ, 2099 MatFDColoringApply_SeqAIJ}; 2100 2101 EXTERN int MatUseSuperLU_SeqAIJ(Mat); 2102 EXTERN int MatUseEssl_SeqAIJ(Mat); 2103 EXTERN int MatUseLUSOL_SeqAIJ(Mat); 2104 EXTERN int MatUseMatlab_SeqAIJ(Mat); 2105 EXTERN int MatUseDXML_SeqAIJ(Mat); 2106 2107 EXTERN_C_BEGIN 2108 #undef __FUNCT__ 2109 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2110 2111 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 2112 { 2113 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2114 int i,nz,n; 2115 2116 PetscFunctionBegin; 2117 2118 nz = aij->maxnz; 2119 n = mat->n; 2120 for (i=0; i<nz; i++) { 2121 aij->j[i] = indices[i]; 2122 } 2123 aij->nz = nz; 2124 for (i=0; i<n; i++) { 2125 aij->ilen[i] = aij->imax[i]; 2126 } 2127 2128 PetscFunctionReturn(0); 2129 } 2130 EXTERN_C_END 2131 2132 #undef __FUNCT__ 2133 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2134 /*@ 2135 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2136 in the matrix. 2137 2138 Input Parameters: 2139 + mat - the SeqAIJ matrix 2140 - indices - the column indices 2141 2142 Level: advanced 2143 2144 Notes: 2145 This can be called if you have precomputed the nonzero structure of the 2146 matrix and want to provide it to the matrix object to improve the performance 2147 of the MatSetValues() operation. 2148 2149 You MUST have set the correct numbers of nonzeros per row in the call to 2150 MatCreateSeqAIJ(). 2151 2152 MUST be called before any calls to MatSetValues(); 2153 2154 The indices should start with zero, not one. 2155 2156 @*/ 2157 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 2158 { 2159 int ierr,(*f)(Mat,int *); 2160 2161 PetscFunctionBegin; 2162 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2163 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)())&f);CHKERRQ(ierr); 2164 if (f) { 2165 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2166 } else { 2167 SETERRQ(1,"Wrong type of matrix to set column indices"); 2168 } 2169 PetscFunctionReturn(0); 2170 } 2171 2172 /* ----------------------------------------------------------------------------------------*/ 2173 2174 EXTERN_C_BEGIN 2175 #undef __FUNCT__ 2176 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2177 int MatStoreValues_SeqAIJ(Mat mat) 2178 { 2179 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2180 int nz = aij->i[mat->m]+aij->indexshift,ierr; 2181 2182 PetscFunctionBegin; 2183 if (aij->nonew != 1) { 2184 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2185 } 2186 2187 /* allocate space for values if not already there */ 2188 if (!aij->saved_values) { 2189 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2190 } 2191 2192 /* copy values over */ 2193 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2194 PetscFunctionReturn(0); 2195 } 2196 EXTERN_C_END 2197 2198 #undef __FUNCT__ 2199 #define __FUNCT__ "MatStoreValues" 2200 /*@ 2201 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2202 example, reuse of the linear part of a Jacobian, while recomputing the 2203 nonlinear portion. 2204 2205 Collect on Mat 2206 2207 Input Parameters: 2208 . mat - the matrix (currently on AIJ matrices support this option) 2209 2210 Level: advanced 2211 2212 Common Usage, with SNESSolve(): 2213 $ Create Jacobian matrix 2214 $ Set linear terms into matrix 2215 $ Apply boundary conditions to matrix, at this time matrix must have 2216 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2217 $ boundary conditions again will not change the nonzero structure 2218 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2219 $ ierr = MatStoreValues(mat); 2220 $ Call SNESSetJacobian() with matrix 2221 $ In your Jacobian routine 2222 $ ierr = MatRetrieveValues(mat); 2223 $ Set nonlinear terms in matrix 2224 2225 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2226 $ // build linear portion of Jacobian 2227 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2228 $ ierr = MatStoreValues(mat); 2229 $ loop over nonlinear iterations 2230 $ ierr = MatRetrieveValues(mat); 2231 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2232 $ // call MatAssemblyBegin/End() on matrix 2233 $ Solve linear system with Jacobian 2234 $ endloop 2235 2236 Notes: 2237 Matrix must already be assemblied before calling this routine 2238 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2239 calling this routine. 2240 2241 .seealso: MatRetrieveValues() 2242 2243 @*/ 2244 int MatStoreValues(Mat mat) 2245 { 2246 int ierr,(*f)(Mat); 2247 2248 PetscFunctionBegin; 2249 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2250 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2251 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2252 2253 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)())&f);CHKERRQ(ierr); 2254 if (f) { 2255 ierr = (*f)(mat);CHKERRQ(ierr); 2256 } else { 2257 SETERRQ(1,"Wrong type of matrix to store values"); 2258 } 2259 PetscFunctionReturn(0); 2260 } 2261 2262 EXTERN_C_BEGIN 2263 #undef __FUNCT__ 2264 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2265 int MatRetrieveValues_SeqAIJ(Mat mat) 2266 { 2267 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2268 int nz = aij->i[mat->m]+aij->indexshift,ierr; 2269 2270 PetscFunctionBegin; 2271 if (aij->nonew != 1) { 2272 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2273 } 2274 if (!aij->saved_values) { 2275 SETERRQ(1,"Must call MatStoreValues(A);first"); 2276 } 2277 2278 /* copy values over */ 2279 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2280 PetscFunctionReturn(0); 2281 } 2282 EXTERN_C_END 2283 2284 #undef __FUNCT__ 2285 #define __FUNCT__ "MatRetrieveValues" 2286 /*@ 2287 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2288 example, reuse of the linear part of a Jacobian, while recomputing the 2289 nonlinear portion. 2290 2291 Collect on Mat 2292 2293 Input Parameters: 2294 . mat - the matrix (currently on AIJ matrices support this option) 2295 2296 Level: advanced 2297 2298 .seealso: MatStoreValues() 2299 2300 @*/ 2301 int MatRetrieveValues(Mat mat) 2302 { 2303 int ierr,(*f)(Mat); 2304 2305 PetscFunctionBegin; 2306 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2307 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2308 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2309 2310 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)())&f);CHKERRQ(ierr); 2311 if (f) { 2312 ierr = (*f)(mat);CHKERRQ(ierr); 2313 } else { 2314 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2315 } 2316 PetscFunctionReturn(0); 2317 } 2318 2319 /* 2320 This allows SeqAIJ matrices to be passed to the matlab engine 2321 */ 2322 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2323 #include "engine.h" /* Matlab include file */ 2324 #include "mex.h" /* Matlab include file */ 2325 EXTERN_C_BEGIN 2326 #undef __FUNCT__ 2327 #define __FUNCT__ "MatMatlabEnginePut_SeqAIJ" 2328 int MatMatlabEnginePut_SeqAIJ(PetscObject obj,void *engine) 2329 { 2330 int ierr,i,*ai,*aj; 2331 Mat B = (Mat)obj; 2332 PetscScalar *array; 2333 mxArray *mat; 2334 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2335 2336 PetscFunctionBegin; 2337 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2338 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2339 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2340 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2341 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2342 2343 /* Matlab indices start at 0 for sparse (what a surprise) */ 2344 if (aij->indexshift) { 2345 for (i=0; i<B->m+1; i++) { 2346 ai[i]--; 2347 } 2348 for (i=0; i<aij->nz; i++) { 2349 aj[i]--; 2350 } 2351 } 2352 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2353 mxSetName(mat,obj->name); 2354 engPutArray((Engine *)engine,mat); 2355 PetscFunctionReturn(0); 2356 } 2357 EXTERN_C_END 2358 2359 EXTERN_C_BEGIN 2360 #undef __FUNCT__ 2361 #define __FUNCT__ "MatMatlabEngineGet_SeqAIJ" 2362 int MatMatlabEngineGet_SeqAIJ(PetscObject obj,void *engine) 2363 { 2364 int ierr,ii; 2365 Mat mat = (Mat)obj; 2366 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2367 mxArray *mmat; 2368 2369 PetscFunctionBegin; 2370 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2371 aij->indexshift = 0; 2372 2373 mmat = engGetArray((Engine *)engine,obj->name); 2374 2375 aij->nz = (mxGetJc(mmat))[mat->m]; 2376 ierr = PetscMalloc(aij->nz*(sizeof(int)+sizeof(PetscScalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 2377 aij->j = (int*)(aij->a + aij->nz); 2378 aij->i = aij->j + aij->nz; 2379 aij->singlemalloc = PETSC_TRUE; 2380 aij->freedata = PETSC_TRUE; 2381 2382 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2383 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2384 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 2385 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 2386 2387 for (ii=0; ii<mat->m; ii++) { 2388 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 2389 } 2390 2391 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2392 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2393 2394 PetscFunctionReturn(0); 2395 } 2396 EXTERN_C_END 2397 #endif 2398 2399 /* --------------------------------------------------------------------------------*/ 2400 2401 #undef __FUNCT__ 2402 #define __FUNCT__ "MatCreateSeqAIJ" 2403 /*@C 2404 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2405 (the default parallel PETSc format). For good matrix assembly performance 2406 the user should preallocate the matrix storage by setting the parameter nz 2407 (or the array nnz). By setting these parameters accurately, performance 2408 during matrix assembly can be increased by more than a factor of 50. 2409 2410 Collective on MPI_Comm 2411 2412 Input Parameters: 2413 + comm - MPI communicator, set to PETSC_COMM_SELF 2414 . m - number of rows 2415 . n - number of columns 2416 . nz - number of nonzeros per row (same for all rows) 2417 - nnz - array containing the number of nonzeros in the various rows 2418 (possibly different for each row) or PETSC_NULL 2419 2420 Output Parameter: 2421 . A - the matrix 2422 2423 Notes: 2424 The AIJ format (also called the Yale sparse matrix format or 2425 compressed row storage), is fully compatible with standard Fortran 77 2426 storage. That is, the stored row and column indices can begin at 2427 either one (as in Fortran) or zero. See the users' manual for details. 2428 2429 Specify the preallocated storage with either nz or nnz (not both). 2430 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2431 allocation. For large problems you MUST preallocate memory or you 2432 will get TERRIBLE performance, see the users' manual chapter on matrices. 2433 2434 By default, this format uses inodes (identical nodes) when possible, to 2435 improve numerical efficiency of matrix-vector products and solves. We 2436 search for consecutive rows with the same nonzero structure, thereby 2437 reusing matrix information to achieve increased efficiency. 2438 2439 Options Database Keys: 2440 + -mat_aij_no_inode - Do not use inodes 2441 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2442 - -mat_aij_oneindex - Internally use indexing starting at 1 2443 rather than 0. Note that when calling MatSetValues(), 2444 the user still MUST index entries starting at 0! 2445 2446 Level: intermediate 2447 2448 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2449 2450 @*/ 2451 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz,Mat *A) 2452 { 2453 int ierr; 2454 2455 PetscFunctionBegin; 2456 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2457 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2458 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2459 PetscFunctionReturn(0); 2460 } 2461 2462 #undef __FUNCT__ 2463 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2464 /*@C 2465 MatSeqAIJSetPreallocation - For good matrix assembly performance 2466 the user should preallocate the matrix storage by setting the parameter nz 2467 (or the array nnz). By setting these parameters accurately, performance 2468 during matrix assembly can be increased by more than a factor of 50. 2469 2470 Collective on MPI_Comm 2471 2472 Input Parameters: 2473 + comm - MPI communicator, set to PETSC_COMM_SELF 2474 . m - number of rows 2475 . n - number of columns 2476 . nz - number of nonzeros per row (same for all rows) 2477 - nnz - array containing the number of nonzeros in the various rows 2478 (possibly different for each row) or PETSC_NULL 2479 2480 Output Parameter: 2481 . A - the matrix 2482 2483 Notes: 2484 The AIJ format (also called the Yale sparse matrix format or 2485 compressed row storage), is fully compatible with standard Fortran 77 2486 storage. That is, the stored row and column indices can begin at 2487 either one (as in Fortran) or zero. See the users' manual for details. 2488 2489 Specify the preallocated storage with either nz or nnz (not both). 2490 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2491 allocation. For large problems you MUST preallocate memory or you 2492 will get TERRIBLE performance, see the users' manual chapter on matrices. 2493 2494 By default, this format uses inodes (identical nodes) when possible, to 2495 improve numerical efficiency of matrix-vector products and solves. We 2496 search for consecutive rows with the same nonzero structure, thereby 2497 reusing matrix information to achieve increased efficiency. 2498 2499 Options Database Keys: 2500 + -mat_aij_no_inode - Do not use inodes 2501 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2502 - -mat_aij_oneindex - Internally use indexing starting at 1 2503 rather than 0. Note that when calling MatSetValues(), 2504 the user still MUST index entries starting at 0! 2505 2506 Level: intermediate 2507 2508 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2509 2510 @*/ 2511 int MatSeqAIJSetPreallocation(Mat B,int nz,int *nnz) 2512 { 2513 Mat_SeqAIJ *b; 2514 int i,len,ierr; 2515 PetscTruth flg2; 2516 2517 PetscFunctionBegin; 2518 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flg2);CHKERRQ(ierr); 2519 if (!flg2) PetscFunctionReturn(0); 2520 2521 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2522 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2523 if (nnz) { 2524 for (i=0; i<B->m; i++) { 2525 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2526 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); 2527 } 2528 } 2529 2530 B->preallocated = PETSC_TRUE; 2531 b = (Mat_SeqAIJ*)B->data; 2532 2533 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2534 if (!nnz) { 2535 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2536 else if (nz <= 0) nz = 1; 2537 for (i=0; i<B->m; i++) b->imax[i] = nz; 2538 nz = nz*B->m; 2539 } else { 2540 nz = 0; 2541 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2542 } 2543 2544 /* allocate the matrix space */ 2545 len = nz*(sizeof(int) + sizeof(PetscScalar)) + (B->m+1)*sizeof(int); 2546 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2547 b->j = (int*)(b->a + nz); 2548 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2549 b->i = b->j + nz; 2550 b->singlemalloc = PETSC_TRUE; 2551 b->freedata = PETSC_TRUE; 2552 2553 b->i[0] = -b->indexshift; 2554 for (i=1; i<B->m+1; i++) { 2555 b->i[i] = b->i[i-1] + b->imax[i-1]; 2556 } 2557 2558 /* b->ilen will count nonzeros in each row so far. */ 2559 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2560 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2561 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2562 2563 b->nz = 0; 2564 b->maxnz = nz; 2565 B->info.nz_unneeded = (double)b->maxnz; 2566 PetscFunctionReturn(0); 2567 } 2568 2569 EXTERN_C_BEGIN 2570 #undef __FUNCT__ 2571 #define __FUNCT__ "MatCreate_SeqAIJ" 2572 int MatCreate_SeqAIJ(Mat B) 2573 { 2574 Mat_SeqAIJ *b; 2575 int ierr,size; 2576 PetscTruth flg; 2577 2578 PetscFunctionBegin; 2579 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2580 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2581 2582 B->m = B->M = PetscMax(B->m,B->M); 2583 B->n = B->N = PetscMax(B->n,B->N); 2584 2585 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2586 B->data = (void*)b; 2587 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2588 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2589 B->factor = 0; 2590 B->lupivotthreshold = 1.0; 2591 B->mapping = 0; 2592 ierr = PetscOptionsGetReal(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2593 ierr = PetscOptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2594 b->row = 0; 2595 b->col = 0; 2596 b->icol = 0; 2597 b->indexshift = 0; 2598 b->reallocs = 0; 2599 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_oneindex",&flg);CHKERRQ(ierr); 2600 if (flg) b->indexshift = -1; 2601 2602 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2603 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2604 2605 b->sorted = PETSC_FALSE; 2606 b->ignorezeroentries = PETSC_FALSE; 2607 b->roworiented = PETSC_TRUE; 2608 b->nonew = 0; 2609 b->diag = 0; 2610 b->solve_work = 0; 2611 b->spptr = 0; 2612 b->inode.use = PETSC_TRUE; 2613 b->inode.node_count = 0; 2614 b->inode.size = 0; 2615 b->inode.limit = 5; 2616 b->inode.max_limit = 5; 2617 b->saved_values = 0; 2618 b->idiag = 0; 2619 b->ssor = 0; 2620 b->keepzeroedrows = PETSC_FALSE; 2621 2622 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2623 2624 #if defined(PETSC_HAVE_SUPERLU) 2625 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_superlu",&flg);CHKERRQ(ierr); 2626 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B);CHKERRQ(ierr); } 2627 #endif 2628 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_essl",&flg);CHKERRQ(ierr); 2629 if (flg) { ierr = MatUseEssl_SeqAIJ(B);CHKERRQ(ierr); } 2630 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_lusol",&flg);CHKERRQ(ierr); 2631 if (flg) { ierr = MatUseLUSOL_SeqAIJ(B);CHKERRQ(ierr); } 2632 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_matlab",&flg);CHKERRQ(ierr); 2633 if (flg) {ierr = MatUseMatlab_SeqAIJ(B);CHKERRQ(ierr);} 2634 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_dxml",&flg);CHKERRQ(ierr); 2635 if (flg) { 2636 if (!b->indexshift) SETERRQ(PETSC_ERR_LIB,"need -mat_aij_oneindex with -mat_aij_dxml"); 2637 ierr = MatUseDXML_SeqAIJ(B);CHKERRQ(ierr); 2638 } 2639 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2640 "MatSeqAIJSetColumnIndices_SeqAIJ", 2641 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2642 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2643 "MatStoreValues_SeqAIJ", 2644 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2645 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2646 "MatRetrieveValues_SeqAIJ", 2647 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2648 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2649 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2650 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2651 #endif 2652 PetscFunctionReturn(0); 2653 } 2654 EXTERN_C_END 2655 2656 #undef __FUNCT__ 2657 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2658 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2659 { 2660 Mat C; 2661 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2662 int i,len,m = A->m,shift = a->indexshift,ierr; 2663 2664 PetscFunctionBegin; 2665 *B = 0; 2666 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2667 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 2668 c = (Mat_SeqAIJ*)C->data; 2669 2670 C->factor = A->factor; 2671 c->row = 0; 2672 c->col = 0; 2673 c->icol = 0; 2674 c->indexshift = shift; 2675 c->keepzeroedrows = a->keepzeroedrows; 2676 C->assembled = PETSC_TRUE; 2677 2678 C->M = A->m; 2679 C->N = A->n; 2680 2681 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2682 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2683 for (i=0; i<m; i++) { 2684 c->imax[i] = a->imax[i]; 2685 c->ilen[i] = a->ilen[i]; 2686 } 2687 2688 /* allocate the matrix space */ 2689 c->singlemalloc = PETSC_TRUE; 2690 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(PetscScalar)+sizeof(int)); 2691 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2692 c->j = (int*)(c->a + a->i[m] + shift); 2693 c->i = c->j + a->i[m] + shift; 2694 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2695 if (m > 0) { 2696 ierr = PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int));CHKERRQ(ierr); 2697 if (cpvalues == MAT_COPY_VALUES) { 2698 ierr = PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(PetscScalar));CHKERRQ(ierr); 2699 } else { 2700 ierr = PetscMemzero(c->a,(a->i[m]+shift)*sizeof(PetscScalar));CHKERRQ(ierr); 2701 } 2702 } 2703 2704 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2705 c->sorted = a->sorted; 2706 c->roworiented = a->roworiented; 2707 c->nonew = a->nonew; 2708 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2709 c->saved_values = 0; 2710 c->idiag = 0; 2711 c->ssor = 0; 2712 c->ignorezeroentries = a->ignorezeroentries; 2713 c->freedata = PETSC_TRUE; 2714 2715 if (a->diag) { 2716 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2717 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2718 for (i=0; i<m; i++) { 2719 c->diag[i] = a->diag[i]; 2720 } 2721 } else c->diag = 0; 2722 c->inode.use = a->inode.use; 2723 c->inode.limit = a->inode.limit; 2724 c->inode.max_limit = a->inode.max_limit; 2725 if (a->inode.size){ 2726 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2727 c->inode.node_count = a->inode.node_count; 2728 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2729 } else { 2730 c->inode.size = 0; 2731 c->inode.node_count = 0; 2732 } 2733 c->nz = a->nz; 2734 c->maxnz = a->maxnz; 2735 c->solve_work = 0; 2736 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2737 C->preallocated = PETSC_TRUE; 2738 2739 *B = C; 2740 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2741 PetscFunctionReturn(0); 2742 } 2743 2744 EXTERN_C_BEGIN 2745 #undef __FUNCT__ 2746 #define __FUNCT__ "MatLoad_SeqAIJ" 2747 int MatLoad_SeqAIJ(PetscViewer viewer,MatType type,Mat *A) 2748 { 2749 Mat_SeqAIJ *a; 2750 Mat B; 2751 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N,shift; 2752 MPI_Comm comm; 2753 2754 PetscFunctionBegin; 2755 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2756 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2757 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2758 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2759 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2760 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2761 M = header[1]; N = header[2]; nz = header[3]; 2762 2763 if (nz < 0) { 2764 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2765 } 2766 2767 /* read in row lengths */ 2768 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2769 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2770 2771 /* create our matrix */ 2772 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A);CHKERRQ(ierr); 2773 B = *A; 2774 a = (Mat_SeqAIJ*)B->data; 2775 shift = a->indexshift; 2776 2777 /* read in column indices and adjust for Fortran indexing*/ 2778 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2779 if (shift) { 2780 for (i=0; i<nz; i++) { 2781 a->j[i] += 1; 2782 } 2783 } 2784 2785 /* read in nonzero values */ 2786 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2787 2788 /* set matrix "i" values */ 2789 a->i[0] = -shift; 2790 for (i=1; i<= M; i++) { 2791 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2792 a->ilen[i-1] = rowlengths[i-1]; 2793 } 2794 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2795 2796 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2797 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2798 PetscFunctionReturn(0); 2799 } 2800 EXTERN_C_END 2801 2802 #undef __FUNCT__ 2803 #define __FUNCT__ "MatEqual_SeqAIJ" 2804 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2805 { 2806 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2807 int ierr; 2808 PetscTruth flag; 2809 2810 PetscFunctionBegin; 2811 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flag);CHKERRQ(ierr); 2812 if (!flag) SETERRQ(PETSC_ERR_ARG_INCOMP,"Matrices must be same type"); 2813 2814 /* If the matrix dimensions are not equal,or no of nonzeros or shift */ 2815 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)|| (a->indexshift != b->indexshift)) { 2816 *flg = PETSC_FALSE; 2817 PetscFunctionReturn(0); 2818 } 2819 2820 /* if the a->i are the same */ 2821 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2822 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2823 2824 /* if a->j are the same */ 2825 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2826 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2827 2828 /* if a->a are the same */ 2829 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2830 2831 PetscFunctionReturn(0); 2832 2833 } 2834 2835 #undef __FUNCT__ 2836 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2837 /*@C 2838 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2839 provided by the user. 2840 2841 Coolective on MPI_Comm 2842 2843 Input Parameters: 2844 + comm - must be an MPI communicator of size 1 2845 . m - number of rows 2846 . n - number of columns 2847 . i - row indices 2848 . j - column indices 2849 - a - matrix values 2850 2851 Output Parameter: 2852 . mat - the matrix 2853 2854 Level: intermediate 2855 2856 Notes: 2857 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2858 once the matrix is destroyed 2859 2860 You cannot set new nonzero locations into this matrix, that will generate an error. 2861 2862 The i and j indices can be either 0- or 1 based 2863 2864 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2865 2866 @*/ 2867 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,PetscScalar *a,Mat *mat) 2868 { 2869 int ierr,ii; 2870 Mat_SeqAIJ *aij; 2871 2872 PetscFunctionBegin; 2873 ierr = MatCreateSeqAIJ(comm,m,n,0,0,mat);CHKERRQ(ierr); 2874 aij = (Mat_SeqAIJ*)(*mat)->data; 2875 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2876 2877 if (i[0] == 1) { 2878 aij->indexshift = -1; 2879 } else if (i[0]) { 2880 SETERRQ(1,"i (row indices) do not start with 0 or 1"); 2881 } 2882 aij->i = i; 2883 aij->j = j; 2884 aij->a = a; 2885 aij->singlemalloc = PETSC_FALSE; 2886 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2887 aij->freedata = PETSC_FALSE; 2888 2889 for (ii=0; ii<m; ii++) { 2890 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2891 #if defined(PETSC_USE_BOPT_g) 2892 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]); 2893 #endif 2894 } 2895 #if defined(PETSC_USE_BOPT_g) 2896 for (ii=0; ii<aij->i[m]; ii++) { 2897 if (j[ii] < -aij->indexshift) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 2898 if (j[ii] > n - 1 -aij->indexshift) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 2899 } 2900 #endif 2901 2902 /* changes indices to start at 0 */ 2903 if (i[0]) { 2904 aij->indexshift = 0; 2905 for (ii=0; ii<m; ii++) { 2906 i[ii]--; 2907 } 2908 for (ii=0; ii<i[m]; ii++) { 2909 j[ii]--; 2910 } 2911 } 2912 2913 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2914 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2915 PetscFunctionReturn(0); 2916 } 2917 2918 #undef __FUNCT__ 2919 #define __FUNCT__ "MatSetColoring_SeqAIJ" 2920 int MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 2921 { 2922 int ierr; 2923 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2924 2925 PetscFunctionBegin; 2926 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 2927 a->coloring = coloring; 2928 PetscFunctionReturn(0); 2929 } 2930 2931 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2932 EXTERN_C_BEGIN 2933 #include "adic_utils.h" 2934 EXTERN_C_END 2935 2936 #undef __FUNCT__ 2937 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 2938 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 2939 { 2940 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2941 int m = A->m,*ii = a->i,*jj = a->j,nz,i,*color,j,nlen; 2942 PetscScalar *v = a->a,*values; 2943 char *cadvalues = (char *)advalues; 2944 2945 PetscFunctionBegin; 2946 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 2947 nlen = my_AD_GetDerivTypeSize(); 2948 color = a->coloring->colors; 2949 /* loop over rows */ 2950 for (i=0; i<m; i++) { 2951 nz = ii[i+1] - ii[i]; 2952 /* loop over columns putting computed value into matrix */ 2953 values = my_AD_GetGradArray(cadvalues); 2954 for (j=0; j<nz; j++) { 2955 *v++ = values[color[*jj++]]; 2956 } 2957 cadvalues += nlen; /* jump to next row of derivatives */ 2958 } 2959 PetscFunctionReturn(0); 2960 } 2961 2962 #else 2963 2964 #undef __FUNCT__ 2965 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 2966 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 2967 { 2968 PetscFunctionBegin; 2969 SETERRQ(1,"PETSc installed without ADIC"); 2970 } 2971 2972 #endif 2973 2974 #undef __FUNCT__ 2975 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 2976 int MatSetValuesAdifor_SeqAIJ(Mat A,int nl,void *advalues) 2977 { 2978 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2979 int m = A->m,*ii = a->i,*jj = a->j,nz,i,*color,j; 2980 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 2981 2982 PetscFunctionBegin; 2983 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 2984 color = a->coloring->colors; 2985 /* loop over rows */ 2986 for (i=0; i<m; i++) { 2987 nz = ii[i+1] - ii[i]; 2988 /* loop over columns putting computed value into matrix */ 2989 for (j=0; j<nz; j++) { 2990 *v++ = values[color[*jj++]]; 2991 } 2992 values += nl; /* jump to next row of derivatives */ 2993 } 2994 PetscFunctionReturn(0); 2995 } 2996 2997