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