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