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