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