Actual source code: aij.c
petsc-3.4.5 2014-06-29
2: /*
3: Defines the basic matrix operations for the AIJ (compressed row)
4: matrix storage format.
5: */
8: #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/
9: #include <petscblaslapack.h>
10: #include <petscbt.h>
11: #include <petsc-private/kernels/blocktranspose.h>
12: #if defined(PETSC_THREADCOMM_ACTIVE)
13: #include <petscthreadcomm.h>
14: #endif
18: PetscErrorCode MatGetColumnNorms_SeqAIJ(Mat A,NormType type,PetscReal *norms)
19: {
21: PetscInt i,m,n;
22: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data;
25: MatGetSize(A,&m,&n);
26: PetscMemzero(norms,n*sizeof(PetscReal));
27: if (type == NORM_2) {
28: for (i=0; i<aij->i[m]; i++) {
29: norms[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]);
30: }
31: } else if (type == NORM_1) {
32: for (i=0; i<aij->i[m]; i++) {
33: norms[aij->j[i]] += PetscAbsScalar(aij->a[i]);
34: }
35: } else if (type == NORM_INFINITY) {
36: for (i=0; i<aij->i[m]; i++) {
37: norms[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),norms[aij->j[i]]);
38: }
39: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown NormType");
41: if (type == NORM_2) {
42: for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
43: }
44: return(0);
45: }
49: PetscErrorCode MatFindZeroDiagonals_SeqAIJ_Private(Mat A,PetscInt *nrows,PetscInt **zrows)
50: {
51: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
52: const MatScalar *aa = a->a;
53: PetscInt i,m=A->rmap->n,cnt = 0;
54: const PetscInt *jj = a->j,*diag;
55: PetscInt *rows;
56: PetscErrorCode ierr;
59: MatMarkDiagonal_SeqAIJ(A);
60: diag = a->diag;
61: for (i=0; i<m; i++) {
62: if ((jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) {
63: cnt++;
64: }
65: }
66: PetscMalloc(cnt*sizeof(PetscInt),&rows);
67: cnt = 0;
68: for (i=0; i<m; i++) {
69: if ((jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) {
70: rows[cnt++] = i;
71: }
72: }
73: *nrows = cnt;
74: *zrows = rows;
75: return(0);
76: }
80: PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A,IS *zrows)
81: {
82: PetscInt nrows,*rows;
86: *zrows = NULL;
87: MatFindZeroDiagonals_SeqAIJ_Private(A,&nrows,&rows);
88: ISCreateGeneral(PetscObjectComm((PetscObject)A),nrows,rows,PETSC_OWN_POINTER,zrows);
89: return(0);
90: }
94: PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A,IS *keptrows)
95: {
96: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
97: const MatScalar *aa;
98: PetscInt m=A->rmap->n,cnt = 0;
99: const PetscInt *ii;
100: PetscInt n,i,j,*rows;
101: PetscErrorCode ierr;
104: *keptrows = 0;
105: ii = a->i;
106: for (i=0; i<m; i++) {
107: n = ii[i+1] - ii[i];
108: if (!n) {
109: cnt++;
110: goto ok1;
111: }
112: aa = a->a + ii[i];
113: for (j=0; j<n; j++) {
114: if (aa[j] != 0.0) goto ok1;
115: }
116: cnt++;
117: ok1:;
118: }
119: if (!cnt) return(0);
120: PetscMalloc((A->rmap->n-cnt)*sizeof(PetscInt),&rows);
121: cnt = 0;
122: for (i=0; i<m; i++) {
123: n = ii[i+1] - ii[i];
124: if (!n) continue;
125: aa = a->a + ii[i];
126: for (j=0; j<n; j++) {
127: if (aa[j] != 0.0) {
128: rows[cnt++] = i;
129: break;
130: }
131: }
132: }
133: ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,keptrows);
134: return(0);
135: }
139: PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is)
140: {
142: Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data;
143: PetscInt i,*diag, m = Y->rmap->n;
144: MatScalar *aa = aij->a;
145: PetscScalar *v;
146: PetscBool missing;
149: if (Y->assembled) {
150: MatMissingDiagonal_SeqAIJ(Y,&missing,NULL);
151: if (!missing) {
152: diag = aij->diag;
153: VecGetArray(D,&v);
154: if (is == INSERT_VALUES) {
155: for (i=0; i<m; i++) {
156: aa[diag[i]] = v[i];
157: }
158: } else {
159: for (i=0; i<m; i++) {
160: aa[diag[i]] += v[i];
161: }
162: }
163: VecRestoreArray(D,&v);
164: return(0);
165: }
166: MatSeqAIJInvalidateDiagonal(Y);
167: }
168: MatDiagonalSet_Default(Y,D,is);
169: return(0);
170: }
174: PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *m,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
175: {
176: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
178: PetscInt i,ishift;
181: *m = A->rmap->n;
182: if (!ia) return(0);
183: ishift = 0;
184: if (symmetric && !A->structurally_symmetric) {
185: MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,ishift,oshift,(PetscInt**)ia,(PetscInt**)ja);
186: } else if (oshift == 1) {
187: PetscInt *tia;
188: PetscInt nz = a->i[A->rmap->n];
189: /* malloc space and add 1 to i and j indices */
190: PetscMalloc((A->rmap->n+1)*sizeof(PetscInt),&tia);
191: for (i=0; i<A->rmap->n+1; i++) tia[i] = a->i[i] + 1;
192: *ia = tia;
193: if (ja) {
194: PetscInt *tja;
195: PetscMalloc((nz+1)*sizeof(PetscInt),&tja);
196: for (i=0; i<nz; i++) tja[i] = a->j[i] + 1;
197: *ja = tja;
198: }
199: } else {
200: *ia = a->i;
201: if (ja) *ja = a->j;
202: }
203: return(0);
204: }
208: PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
209: {
213: if (!ia) return(0);
214: if ((symmetric && !A->structurally_symmetric) || oshift == 1) {
215: PetscFree(*ia);
216: if (ja) {PetscFree(*ja);}
217: }
218: return(0);
219: }
223: PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
224: {
225: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
227: PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n;
228: PetscInt nz = a->i[m],row,*jj,mr,col;
231: *nn = n;
232: if (!ia) return(0);
233: if (symmetric) {
234: MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,0,oshift,(PetscInt**)ia,(PetscInt**)ja);
235: } else {
236: PetscMalloc((n+1)*sizeof(PetscInt),&collengths);
237: PetscMemzero(collengths,n*sizeof(PetscInt));
238: PetscMalloc((n+1)*sizeof(PetscInt),&cia);
239: PetscMalloc((nz+1)*sizeof(PetscInt),&cja);
240: jj = a->j;
241: for (i=0; i<nz; i++) {
242: collengths[jj[i]]++;
243: }
244: cia[0] = oshift;
245: for (i=0; i<n; i++) {
246: cia[i+1] = cia[i] + collengths[i];
247: }
248: PetscMemzero(collengths,n*sizeof(PetscInt));
249: jj = a->j;
250: for (row=0; row<m; row++) {
251: mr = a->i[row+1] - a->i[row];
252: for (i=0; i<mr; i++) {
253: col = *jj++;
255: cja[cia[col] + collengths[col]++ - oshift] = row + oshift;
256: }
257: }
258: PetscFree(collengths);
259: *ia = cia; *ja = cja;
260: }
261: return(0);
262: }
266: PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
267: {
271: if (!ia) return(0);
273: PetscFree(*ia);
274: PetscFree(*ja);
275: return(0);
276: }
280: PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[])
281: {
282: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
283: PetscInt *ai = a->i;
287: PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));
288: return(0);
289: }
293: PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
294: {
295: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
296: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N;
297: PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen;
299: PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1;
300: MatScalar *ap,value,*aa = a->a;
301: PetscBool ignorezeroentries = a->ignorezeroentries;
302: PetscBool roworiented = a->roworiented;
306: for (k=0; k<m; k++) { /* loop over added rows */
307: row = im[k];
308: if (row < 0) continue;
309: #if defined(PETSC_USE_DEBUG)
310: if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1);
311: #endif
312: rp = aj + ai[row]; ap = aa + ai[row];
313: rmax = imax[row]; nrow = ailen[row];
314: low = 0;
315: high = nrow;
316: for (l=0; l<n; l++) { /* loop over added columns */
317: if (in[l] < 0) continue;
318: #if defined(PETSC_USE_DEBUG)
319: if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1);
320: #endif
321: col = in[l];
322: if (v) {
323: if (roworiented) {
324: value = v[l + k*n];
325: } else {
326: value = v[k + l*m];
327: }
328: } else {
329: value = 0.;
330: }
331: if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue;
333: if (col <= lastcol) low = 0;
334: else high = nrow;
335: lastcol = col;
336: while (high-low > 5) {
337: t = (low+high)/2;
338: if (rp[t] > col) high = t;
339: else low = t;
340: }
341: for (i=low; i<high; i++) {
342: if (rp[i] > col) break;
343: if (rp[i] == col) {
344: if (is == ADD_VALUES) ap[i] += value;
345: else ap[i] = value;
346: low = i + 1;
347: goto noinsert;
348: }
349: }
350: if (value == 0.0 && ignorezeroentries) goto noinsert;
351: if (nonew == 1) goto noinsert;
352: if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col);
353: MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
354: N = nrow++ - 1; a->nz++; high++;
355: /* shift up all the later entries in this row */
356: for (ii=N; ii>=i; ii--) {
357: rp[ii+1] = rp[ii];
358: ap[ii+1] = ap[ii];
359: }
360: rp[i] = col;
361: ap[i] = value;
362: low = i + 1;
363: noinsert:;
364: }
365: ailen[row] = nrow;
366: }
367: A->same_nonzero = PETSC_FALSE;
368: return(0);
369: }
374: PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[])
375: {
376: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
377: PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j;
378: PetscInt *ai = a->i,*ailen = a->ilen;
379: MatScalar *ap,*aa = a->a;
382: for (k=0; k<m; k++) { /* loop over rows */
383: row = im[k];
384: if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */
385: if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1);
386: rp = aj + ai[row]; ap = aa + ai[row];
387: nrow = ailen[row];
388: for (l=0; l<n; l++) { /* loop over columns */
389: if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */
390: if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1);
391: col = in[l];
392: high = nrow; low = 0; /* assume unsorted */
393: while (high-low > 5) {
394: t = (low+high)/2;
395: if (rp[t] > col) high = t;
396: else low = t;
397: }
398: for (i=low; i<high; i++) {
399: if (rp[i] > col) break;
400: if (rp[i] == col) {
401: *v++ = ap[i];
402: goto finished;
403: }
404: }
405: *v++ = 0.0;
406: finished:;
407: }
408: }
409: return(0);
410: }
415: PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer)
416: {
417: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
419: PetscInt i,*col_lens;
420: int fd;
421: FILE *file;
424: PetscViewerBinaryGetDescriptor(viewer,&fd);
425: PetscMalloc((4+A->rmap->n)*sizeof(PetscInt),&col_lens);
427: col_lens[0] = MAT_FILE_CLASSID;
428: col_lens[1] = A->rmap->n;
429: col_lens[2] = A->cmap->n;
430: col_lens[3] = a->nz;
432: /* store lengths of each row and write (including header) to file */
433: for (i=0; i<A->rmap->n; i++) {
434: col_lens[4+i] = a->i[i+1] - a->i[i];
435: }
436: PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);
437: PetscFree(col_lens);
439: /* store column indices (zero start index) */
440: PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);
442: /* store nonzero values */
443: PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);
445: PetscViewerBinaryGetInfoPointer(viewer,&file);
446: if (file) {
447: fprintf(file,"-matload_block_size %d\n",(int)A->rmap->bs);
448: }
449: return(0);
450: }
452: extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer);
456: PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer)
457: {
458: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
459: PetscErrorCode ierr;
460: PetscInt i,j,m = A->rmap->n,shift=0;
461: const char *name;
462: PetscViewerFormat format;
465: PetscViewerGetFormat(viewer,&format);
466: if (format == PETSC_VIEWER_ASCII_MATLAB) {
467: PetscInt nofinalvalue = 0;
468: if (m && ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-!shift))) {
469: nofinalvalue = 1;
470: }
471: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
472: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap->n);
473: PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);
474: PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);
475: PetscViewerASCIIPrintf(viewer,"zzz = [\n");
477: for (i=0; i<m; i++) {
478: for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) {
479: #if defined(PETSC_USE_COMPLEX)
480: PetscViewerASCIIPrintf(viewer,"%D %D %18.16e + %18.16ei \n",i+1,a->j[j]+!shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));
481: #else
482: PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);
483: #endif
484: }
485: }
486: if (nofinalvalue) {
487: PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap->n,0.0);
488: }
489: PetscObjectGetName((PetscObject)A,&name);
490: PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);
491: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
492: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) {
493: return(0);
494: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
495: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
496: PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");
497: for (i=0; i<m; i++) {
498: PetscViewerASCIIPrintf(viewer,"row %D:",i);
499: for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) {
500: #if defined(PETSC_USE_COMPLEX)
501: if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) {
502: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));
503: } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) {
504: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));
505: } else if (PetscRealPart(a->a[j]) != 0.0) {
506: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));
507: }
508: #else
509: if (a->a[j] != 0.0) {PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,(double)a->a[j]);}
510: #endif
511: }
512: PetscViewerASCIIPrintf(viewer,"\n");
513: }
514: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
515: } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
516: PetscInt nzd=0,fshift=1,*sptr;
517: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
518: PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");
519: PetscMalloc((m+1)*sizeof(PetscInt),&sptr);
520: for (i=0; i<m; i++) {
521: sptr[i] = nzd+1;
522: for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) {
523: if (a->j[j] >= i) {
524: #if defined(PETSC_USE_COMPLEX)
525: if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++;
526: #else
527: if (a->a[j] != 0.0) nzd++;
528: #endif
529: }
530: }
531: }
532: sptr[m] = nzd+1;
533: PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);
534: for (i=0; i<m+1; i+=6) {
535: if (i+4<m) {
536: PetscViewerASCIIPrintf(viewer," %D %D %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5]);
537: } else if (i+3<m) {
538: PetscViewerASCIIPrintf(viewer," %D %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4]);
539: } else if (i+2<m) {
540: PetscViewerASCIIPrintf(viewer," %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]);
541: } else if (i+1<m) {
542: PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);
543: } else if (i<m) {
544: PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);
545: } else {
546: PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);
547: }
548: }
549: PetscViewerASCIIPrintf(viewer,"\n");
550: PetscFree(sptr);
551: for (i=0; i<m; i++) {
552: for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) {
553: if (a->j[j] >= i) {PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);}
554: }
555: PetscViewerASCIIPrintf(viewer,"\n");
556: }
557: PetscViewerASCIIPrintf(viewer,"\n");
558: for (i=0; i<m; i++) {
559: for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) {
560: if (a->j[j] >= i) {
561: #if defined(PETSC_USE_COMPLEX)
562: if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) {
563: PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));
564: }
565: #else
566: if (a->a[j] != 0.0) {PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);}
567: #endif
568: }
569: }
570: PetscViewerASCIIPrintf(viewer,"\n");
571: }
572: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
573: } else if (format == PETSC_VIEWER_ASCII_DENSE) {
574: PetscInt cnt = 0,jcnt;
575: PetscScalar value;
577: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
578: PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");
579: for (i=0; i<m; i++) {
580: jcnt = 0;
581: for (j=0; j<A->cmap->n; j++) {
582: if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) {
583: value = a->a[cnt++];
584: jcnt++;
585: } else {
586: value = 0.0;
587: }
588: #if defined(PETSC_USE_COMPLEX)
589: PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));
590: #else
591: PetscViewerASCIIPrintf(viewer," %7.5e ",value);
592: #endif
593: }
594: PetscViewerASCIIPrintf(viewer,"\n");
595: }
596: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
597: } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) {
598: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
599: PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");
600: #if defined(PETSC_USE_COMPLEX)
601: PetscViewerASCIIPrintf(viewer,"%%matrix complex general\n");
602: #else
603: PetscViewerASCIIPrintf(viewer,"%%matrix real general\n");
604: #endif
605: PetscViewerASCIIPrintf(viewer,"%D %D %D\n", m, A->cmap->n, a->nz);
606: for (i=0; i<m; i++) {
607: for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) {
608: #if defined(PETSC_USE_COMPLEX)
609: if (PetscImaginaryPart(a->a[j]) > 0.0) {
610: PetscViewerASCIIPrintf(viewer,"%D %D, %g %g\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));
611: } else if (PetscImaginaryPart(a->a[j]) < 0.0) {
612: PetscViewerASCIIPrintf(viewer,"%D %D, %g -%g\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));
613: } else {
614: PetscViewerASCIIPrintf(viewer,"%D %D, %g\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]));
615: }
616: #else
617: PetscViewerASCIIPrintf(viewer,"%D %D %g\n", i+shift, a->j[j]+shift, (double)a->a[j]);
618: #endif
619: }
620: }
621: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
622: } else {
623: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
624: PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");
625: if (A->factortype) {
626: for (i=0; i<m; i++) {
627: PetscViewerASCIIPrintf(viewer,"row %D:",i);
628: /* L part */
629: for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) {
630: #if defined(PETSC_USE_COMPLEX)
631: if (PetscImaginaryPart(a->a[j]) > 0.0) {
632: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));
633: } else if (PetscImaginaryPart(a->a[j]) < 0.0) {
634: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));
635: } else {
636: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));
637: }
638: #else
639: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,(double)a->a[j]);
640: #endif
641: }
642: /* diagonal */
643: j = a->diag[i];
644: #if defined(PETSC_USE_COMPLEX)
645: if (PetscImaginaryPart(a->a[j]) > 0.0) {
646: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(1.0/a->a[j]),PetscImaginaryPart(1.0/a->a[j]));
647: } else if (PetscImaginaryPart(a->a[j]) < 0.0) {
648: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(1.0/a->a[j]),-PetscImaginaryPart(1.0/a->a[j]));
649: } else {
650: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(1.0/a->a[j]));
651: }
652: #else
653: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,(double)1.0/a->a[j]);
654: #endif
656: /* U part */
657: for (j=a->diag[i+1]+1+shift; j<a->diag[i]+shift; j++) {
658: #if defined(PETSC_USE_COMPLEX)
659: if (PetscImaginaryPart(a->a[j]) > 0.0) {
660: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));
661: } else if (PetscImaginaryPart(a->a[j]) < 0.0) {
662: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));
663: } else {
664: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));
665: }
666: #else
667: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,(double)a->a[j]);
668: #endif
669: }
670: PetscViewerASCIIPrintf(viewer,"\n");
671: }
672: } else {
673: for (i=0; i<m; i++) {
674: PetscViewerASCIIPrintf(viewer,"row %D:",i);
675: for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) {
676: #if defined(PETSC_USE_COMPLEX)
677: if (PetscImaginaryPart(a->a[j]) > 0.0) {
678: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));
679: } else if (PetscImaginaryPart(a->a[j]) < 0.0) {
680: PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));
681: } else {
682: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));
683: }
684: #else
685: PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,(double)a->a[j]);
686: #endif
687: }
688: PetscViewerASCIIPrintf(viewer,"\n");
689: }
690: }
691: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
692: }
693: PetscViewerFlush(viewer);
694: return(0);
695: }
697: #include <petscdraw.h>
700: PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa)
701: {
702: Mat A = (Mat) Aa;
703: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
704: PetscErrorCode ierr;
705: PetscInt i,j,m = A->rmap->n,color;
706: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0;
707: PetscViewer viewer;
708: PetscViewerFormat format;
711: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
712: PetscViewerGetFormat(viewer,&format);
714: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
715: /* loop over matrix elements drawing boxes */
717: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
718: /* Blue for negative, Cyan for zero and Red for positive */
719: color = PETSC_DRAW_BLUE;
720: for (i=0; i<m; i++) {
721: y_l = m - i - 1.0; y_r = y_l + 1.0;
722: for (j=a->i[i]; j<a->i[i+1]; j++) {
723: x_l = a->j[j]; x_r = x_l + 1.0;
724: if (PetscRealPart(a->a[j]) >= 0.) continue;
725: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
726: }
727: }
728: color = PETSC_DRAW_CYAN;
729: for (i=0; i<m; i++) {
730: y_l = m - i - 1.0; y_r = y_l + 1.0;
731: for (j=a->i[i]; j<a->i[i+1]; j++) {
732: x_l = a->j[j]; x_r = x_l + 1.0;
733: if (a->a[j] != 0.) continue;
734: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
735: }
736: }
737: color = PETSC_DRAW_RED;
738: for (i=0; i<m; i++) {
739: y_l = m - i - 1.0; y_r = y_l + 1.0;
740: for (j=a->i[i]; j<a->i[i+1]; j++) {
741: x_l = a->j[j]; x_r = x_l + 1.0;
742: if (PetscRealPart(a->a[j]) <= 0.) continue;
743: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
744: }
745: }
746: } else {
747: /* use contour shading to indicate magnitude of values */
748: /* first determine max of all nonzero values */
749: PetscInt nz = a->nz,count;
750: PetscDraw popup;
751: PetscReal scale;
753: for (i=0; i<nz; i++) {
754: if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]);
755: }
756: scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv;
757: PetscDrawGetPopup(draw,&popup);
758: if (popup) {
759: PetscDrawScalePopup(popup,0.0,maxv);
760: }
761: count = 0;
762: for (i=0; i<m; i++) {
763: y_l = m - i - 1.0; y_r = y_l + 1.0;
764: for (j=a->i[i]; j<a->i[i+1]; j++) {
765: x_l = a->j[j]; x_r = x_l + 1.0;
766: color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count]));
767: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
768: count++;
769: }
770: }
771: }
772: return(0);
773: }
775: #include <petscdraw.h>
778: PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer)
779: {
781: PetscDraw draw;
782: PetscReal xr,yr,xl,yl,h,w;
783: PetscBool isnull;
786: PetscViewerDrawGetDraw(viewer,0,&draw);
787: PetscDrawIsNull(draw,&isnull);
788: if (isnull) return(0);
790: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
791: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
792: xr += w; yr += h; xl = -w; yl = -h;
793: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
794: PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);
795: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
796: return(0);
797: }
801: PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer)
802: {
804: PetscBool iascii,isbinary,isdraw;
807: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
808: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
809: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
810: if (iascii) {
811: MatView_SeqAIJ_ASCII(A,viewer);
812: } else if (isbinary) {
813: MatView_SeqAIJ_Binary(A,viewer);
814: } else if (isdraw) {
815: MatView_SeqAIJ_Draw(A,viewer);
816: }
817: MatView_SeqAIJ_Inode(A,viewer);
818: return(0);
819: }
823: PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode)
824: {
825: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
827: PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax;
828: PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0;
829: MatScalar *aa = a->a,*ap;
830: PetscReal ratio = 0.6;
833: if (mode == MAT_FLUSH_ASSEMBLY) return(0);
835: if (m) rmax = ailen[0]; /* determine row with most nonzeros */
836: for (i=1; i<m; i++) {
837: /* move each row back by the amount of empty slots (fshift) before it*/
838: fshift += imax[i-1] - ailen[i-1];
839: rmax = PetscMax(rmax,ailen[i]);
840: if (fshift) {
841: ip = aj + ai[i];
842: ap = aa + ai[i];
843: N = ailen[i];
844: for (j=0; j<N; j++) {
845: ip[j-fshift] = ip[j];
846: ap[j-fshift] = ap[j];
847: }
848: }
849: ai[i] = ai[i-1] + ailen[i-1];
850: }
851: if (m) {
852: fshift += imax[m-1] - ailen[m-1];
853: ai[m] = ai[m-1] + ailen[m-1];
854: }
855: /* reset ilen and imax for each row */
856: for (i=0; i<m; i++) {
857: ailen[i] = imax[i] = ai[i+1] - ai[i];
858: }
859: a->nz = ai[m];
860: if (fshift && a->nounused == -1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D, %D unneeded", m, A->cmap->n, fshift);
862: MatMarkDiagonal_SeqAIJ(A);
863: PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n,fshift,a->nz);
864: PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);
865: PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);
867: A->info.mallocs += a->reallocs;
868: a->reallocs = 0;
869: A->info.nz_unneeded = (double)fshift;
870: a->rmax = rmax;
872: MatCheckCompressedRow(A,&a->compressedrow,a->i,m,ratio);
874: A->same_nonzero = PETSC_TRUE;
876: MatAssemblyEnd_SeqAIJ_Inode(A,mode);
878: MatSeqAIJInvalidateDiagonal(A);
879: return(0);
880: }
884: PetscErrorCode MatRealPart_SeqAIJ(Mat A)
885: {
886: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
887: PetscInt i,nz = a->nz;
888: MatScalar *aa = a->a;
892: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
893: MatSeqAIJInvalidateDiagonal(A);
894: return(0);
895: }
899: PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A)
900: {
901: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
902: PetscInt i,nz = a->nz;
903: MatScalar *aa = a->a;
907: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
908: MatSeqAIJInvalidateDiagonal(A);
909: return(0);
910: }
912: #if defined(PETSC_THREADCOMM_ACTIVE)
913: PetscErrorCode MatZeroEntries_SeqAIJ_Kernel(PetscInt thread_id,Mat A)
914: {
916: PetscInt *trstarts=A->rmap->trstarts;
917: PetscInt n,start,end;
918: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
920: start = trstarts[thread_id];
921: end = trstarts[thread_id+1];
922: n = a->i[end] - a->i[start];
923: PetscMemzero(a->a+a->i[start],n*sizeof(PetscScalar));
924: return 0;
925: }
929: PetscErrorCode MatZeroEntries_SeqAIJ(Mat A)
930: {
934: PetscThreadCommRunKernel(PetscObjectComm((PetscObject)A),(PetscThreadKernel)MatZeroEntries_SeqAIJ_Kernel,1,A);
935: MatSeqAIJInvalidateDiagonal(A);
936: return(0);
937: }
938: #else
941: PetscErrorCode MatZeroEntries_SeqAIJ(Mat A)
942: {
943: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
947: PetscMemzero(a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));
948: MatSeqAIJInvalidateDiagonal(A);
949: return(0);
950: }
951: #endif
955: PetscErrorCode MatDestroy_SeqAIJ(Mat A)
956: {
957: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
961: #if defined(PETSC_USE_LOG)
962: PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap->n,A->cmap->n,a->nz);
963: #endif
964: MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);
965: ISDestroy(&a->row);
966: ISDestroy(&a->col);
967: PetscFree(a->diag);
968: PetscFree(a->ibdiag);
969: PetscFree2(a->imax,a->ilen);
970: PetscFree3(a->idiag,a->mdiag,a->ssor_work);
971: PetscFree(a->solve_work);
972: ISDestroy(&a->icol);
973: PetscFree(a->saved_values);
974: ISColoringDestroy(&a->coloring);
975: PetscFree(a->xtoy);
976: MatDestroy(&a->XtoY);
977: PetscFree2(a->compressedrow.i,a->compressedrow.rindex);
978: PetscFree(a->matmult_abdense);
980: MatDestroy_SeqAIJ_Inode(A);
981: PetscFree(A->data);
983: PetscObjectChangeTypeName((PetscObject)A,0);
984: PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetColumnIndices_C",NULL);
985: PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);
986: PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);
987: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsbaij_C",NULL);
988: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqbaij_C",NULL);
989: PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijperm_C",NULL);
990: PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL);
991: PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL);
992: PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL);
993: PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL);
994: return(0);
995: }
999: PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg)
1000: {
1001: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1005: switch (op) {
1006: case MAT_ROW_ORIENTED:
1007: a->roworiented = flg;
1008: break;
1009: case MAT_KEEP_NONZERO_PATTERN:
1010: a->keepnonzeropattern = flg;
1011: break;
1012: case MAT_NEW_NONZERO_LOCATIONS:
1013: a->nonew = (flg ? 0 : 1);
1014: break;
1015: case MAT_NEW_NONZERO_LOCATION_ERR:
1016: a->nonew = (flg ? -1 : 0);
1017: break;
1018: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1019: a->nonew = (flg ? -2 : 0);
1020: break;
1021: case MAT_UNUSED_NONZERO_LOCATION_ERR:
1022: a->nounused = (flg ? -1 : 0);
1023: break;
1024: case MAT_IGNORE_ZERO_ENTRIES:
1025: a->ignorezeroentries = flg;
1026: break;
1027: case MAT_CHECK_COMPRESSED_ROW:
1028: a->compressedrow.check = flg;
1029: break;
1030: case MAT_SPD:
1031: case MAT_SYMMETRIC:
1032: case MAT_STRUCTURALLY_SYMMETRIC:
1033: case MAT_HERMITIAN:
1034: case MAT_SYMMETRY_ETERNAL:
1035: /* These options are handled directly by MatSetOption() */
1036: break;
1037: case MAT_NEW_DIAGONALS:
1038: case MAT_IGNORE_OFF_PROC_ENTRIES:
1039: case MAT_USE_HASH_TABLE:
1040: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1041: break;
1042: case MAT_USE_INODES:
1043: /* Not an error because MatSetOption_SeqAIJ_Inode handles this one */
1044: break;
1045: default:
1046: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op);
1047: }
1048: MatSetOption_SeqAIJ_Inode(A,op,flg);
1049: return(0);
1050: }
1054: PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v)
1055: {
1056: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1058: PetscInt i,j,n,*ai=a->i,*aj=a->j,nz;
1059: PetscScalar *aa=a->a,*x,zero=0.0;
1062: VecGetLocalSize(v,&n);
1063: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
1065: if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) {
1066: PetscInt *diag=a->diag;
1067: VecGetArray(v,&x);
1068: for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]];
1069: VecRestoreArray(v,&x);
1070: return(0);
1071: }
1073: VecSet(v,zero);
1074: VecGetArray(v,&x);
1075: for (i=0; i<n; i++) {
1076: nz = ai[i+1] - ai[i];
1077: if (!nz) x[i] = 0.0;
1078: for (j=ai[i]; j<ai[i+1]; j++) {
1079: if (aj[j] == i) {
1080: x[i] = aa[j];
1081: break;
1082: }
1083: }
1084: }
1085: VecRestoreArray(v,&x);
1086: return(0);
1087: }
1089: #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h>
1092: PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy)
1093: {
1094: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1095: PetscScalar *x,*y;
1097: PetscInt m = A->rmap->n;
1098: #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ)
1099: MatScalar *v;
1100: PetscScalar alpha;
1101: PetscInt n,i,j,*idx,*ii,*ridx=NULL;
1102: Mat_CompressedRow cprow = a->compressedrow;
1103: PetscBool usecprow = cprow.use;
1104: #endif
1107: if (zz != yy) {VecCopy(zz,yy);}
1108: VecGetArray(xx,&x);
1109: VecGetArray(yy,&y);
1111: #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ)
1112: fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y);
1113: #else
1114: if (usecprow) {
1115: m = cprow.nrows;
1116: ii = cprow.i;
1117: ridx = cprow.rindex;
1118: } else {
1119: ii = a->i;
1120: }
1121: for (i=0; i<m; i++) {
1122: idx = a->j + ii[i];
1123: v = a->a + ii[i];
1124: n = ii[i+1] - ii[i];
1125: if (usecprow) {
1126: alpha = x[ridx[i]];
1127: } else {
1128: alpha = x[i];
1129: }
1130: for (j=0; j<n; j++) y[idx[j]] += alpha*v[j];
1131: }
1132: #endif
1133: PetscLogFlops(2.0*a->nz);
1134: VecRestoreArray(xx,&x);
1135: VecRestoreArray(yy,&y);
1136: return(0);
1137: }
1141: PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy)
1142: {
1146: VecSet(yy,0.0);
1147: MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);
1148: return(0);
1149: }
1151: #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h>
1152: #if defined(PETSC_THREADCOMM_ACTIVE)
1153: PetscErrorCode MatMult_SeqAIJ_Kernel(PetscInt thread_id,Mat A,Vec xx,Vec yy)
1154: {
1155: PetscErrorCode ierr;
1156: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1157: PetscScalar *y;
1158: const PetscScalar *x;
1159: const MatScalar *aa;
1160: PetscInt *trstarts=A->rmap->trstarts;
1161: PetscInt n,start,end,i;
1162: const PetscInt *aj,*ai;
1163: PetscScalar sum;
1165: VecGetArrayRead(xx,&x);
1166: VecGetArray(yy,&y);
1167: start = trstarts[thread_id];
1168: end = trstarts[thread_id+1];
1169: aj = a->j;
1170: aa = a->a;
1171: ai = a->i;
1172: for (i=start; i<end; i++) {
1173: n = ai[i+1] - ai[i];
1174: aj = a->j + ai[i];
1175: aa = a->a + ai[i];
1176: sum = 0.0;
1177: PetscSparseDensePlusDot(sum,x,aa,aj,n);
1178: y[i] = sum;
1179: }
1180: VecRestoreArrayRead(xx,&x);
1181: VecRestoreArray(yy,&y);
1182: return 0;
1183: }
1187: PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy)
1188: {
1189: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1190: PetscScalar *y;
1191: const PetscScalar *x;
1192: const MatScalar *aa;
1193: PetscErrorCode ierr;
1194: PetscInt m=A->rmap->n;
1195: const PetscInt *aj,*ii,*ridx=NULL;
1196: PetscInt n,i,nonzerorow=0;
1197: PetscScalar sum;
1198: PetscBool usecprow=a->compressedrow.use;
1200: #if defined(PETSC_HAVE_PRAGMA_DISJOINT)
1201: #pragma disjoint(*x,*y,*aa)
1202: #endif
1205: aj = a->j;
1206: aa = a->a;
1207: ii = a->i;
1208: if (usecprow) { /* use compressed row format */
1209: VecGetArrayRead(xx,&x);
1210: VecGetArray(yy,&y);
1211: m = a->compressedrow.nrows;
1212: ii = a->compressedrow.i;
1213: ridx = a->compressedrow.rindex;
1214: for (i=0; i<m; i++) {
1215: n = ii[i+1] - ii[i];
1216: aj = a->j + ii[i];
1217: aa = a->a + ii[i];
1218: sum = 0.0;
1219: nonzerorow += (n>0);
1220: PetscSparseDensePlusDot(sum,x,aa,aj,n);
1221: /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */
1222: y[*ridx++] = sum;
1223: }
1224: VecRestoreArrayRead(xx,&x);
1225: VecRestoreArray(yy,&y);
1226: } else { /* do not use compressed row format */
1227: #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ)
1228: fortranmultaij_(&m,x,ii,aj,aa,y);
1229: #else
1230: PetscThreadCommRunKernel(PetscObjectComm((PetscObject)A),(PetscThreadKernel)MatMult_SeqAIJ_Kernel,3,A,xx,yy);
1231: #endif
1232: }
1233: PetscLogFlops(2.0*a->nz - nonzerorow);
1234: return(0);
1235: }
1236: #else
1239: PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy)
1240: {
1241: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1242: PetscScalar *y;
1243: const PetscScalar *x;
1244: const MatScalar *aa;
1245: PetscErrorCode ierr;
1246: PetscInt m=A->rmap->n;
1247: const PetscInt *aj,*ii,*ridx=NULL;
1248: PetscInt n,i,nonzerorow=0;
1249: PetscScalar sum;
1250: PetscBool usecprow=a->compressedrow.use;
1252: #if defined(PETSC_HAVE_PRAGMA_DISJOINT)
1253: #pragma disjoint(*x,*y,*aa)
1254: #endif
1257: VecGetArrayRead(xx,&x);
1258: VecGetArray(yy,&y);
1259: aj = a->j;
1260: aa = a->a;
1261: ii = a->i;
1262: if (usecprow) { /* use compressed row format */
1263: m = a->compressedrow.nrows;
1264: ii = a->compressedrow.i;
1265: ridx = a->compressedrow.rindex;
1266: for (i=0; i<m; i++) {
1267: n = ii[i+1] - ii[i];
1268: aj = a->j + ii[i];
1269: aa = a->a + ii[i];
1270: sum = 0.0;
1271: nonzerorow += (n>0);
1272: PetscSparseDensePlusDot(sum,x,aa,aj,n);
1273: /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */
1274: y[*ridx++] = sum;
1275: }
1276: } else { /* do not use compressed row format */
1277: #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ)
1278: fortranmultaij_(&m,x,ii,aj,aa,y);
1279: #else
1280: #if defined(PETSC_THREADCOMM_ACTIVE)
1281: PetscThreadCommRunKernel(PetscObjectComm((PetscObject)A),(PetscThreadKernel)MatMult_SeqAIJ_Kernel,3,A,xx,yy);
1282: #else
1283: for (i=0; i<m; i++) {
1284: n = ii[i+1] - ii[i];
1285: aj = a->j + ii[i];
1286: aa = a->a + ii[i];
1287: sum = 0.0;
1288: nonzerorow += (n>0);
1289: PetscSparseDensePlusDot(sum,x,aa,aj,n);
1290: y[i] = sum;
1291: }
1292: #endif
1293: #endif
1294: }
1295: PetscLogFlops(2.0*a->nz - nonzerorow);
1296: VecRestoreArrayRead(xx,&x);
1297: VecRestoreArray(yy,&y);
1298: return(0);
1299: }
1300: #endif
1302: #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h>
1305: PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz)
1306: {
1307: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1308: PetscScalar *y,*z;
1309: const PetscScalar *x;
1310: const MatScalar *aa;
1311: PetscErrorCode ierr;
1312: PetscInt m = A->rmap->n,*aj,*ii;
1313: PetscInt n,i,*ridx=NULL;
1314: PetscScalar sum;
1315: PetscBool usecprow=a->compressedrow.use;
1318: VecGetArrayRead(xx,&x);
1319: VecGetArray(yy,&y);
1320: if (zz != yy) {
1321: VecGetArray(zz,&z);
1322: } else {
1323: z = y;
1324: }
1326: aj = a->j;
1327: aa = a->a;
1328: ii = a->i;
1329: if (usecprow) { /* use compressed row format */
1330: if (zz != yy) {
1331: PetscMemcpy(z,y,m*sizeof(PetscScalar));
1332: }
1333: m = a->compressedrow.nrows;
1334: ii = a->compressedrow.i;
1335: ridx = a->compressedrow.rindex;
1336: for (i=0; i<m; i++) {
1337: n = ii[i+1] - ii[i];
1338: aj = a->j + ii[i];
1339: aa = a->a + ii[i];
1340: sum = y[*ridx];
1341: PetscSparseDensePlusDot(sum,x,aa,aj,n);
1342: z[*ridx++] = sum;
1343: }
1344: } else { /* do not use compressed row format */
1345: #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ)
1346: fortranmultaddaij_(&m,x,ii,aj,aa,y,z);
1347: #else
1348: for (i=0; i<m; i++) {
1349: n = ii[i+1] - ii[i];
1350: aj = a->j + ii[i];
1351: aa = a->a + ii[i];
1352: sum = y[i];
1353: PetscSparseDensePlusDot(sum,x,aa,aj,n);
1354: z[i] = sum;
1355: }
1356: #endif
1357: }
1358: PetscLogFlops(2.0*a->nz);
1359: VecRestoreArrayRead(xx,&x);
1360: VecRestoreArray(yy,&y);
1361: if (zz != yy) {
1362: VecRestoreArray(zz,&z);
1363: }
1364: #if defined(PETSC_HAVE_CUSP)
1365: /*
1366: VecView(xx,0);
1367: VecView(zz,0);
1368: MatView(A,0);
1369: */
1370: #endif
1371: return(0);
1372: }
1374: /*
1375: Adds diagonal pointers to sparse matrix structure.
1376: */
1379: PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A)
1380: {
1381: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1383: PetscInt i,j,m = A->rmap->n;
1386: if (!a->diag) {
1387: PetscMalloc(m*sizeof(PetscInt),&a->diag);
1388: PetscLogObjectMemory(A, m*sizeof(PetscInt));
1389: }
1390: for (i=0; i<A->rmap->n; i++) {
1391: a->diag[i] = a->i[i+1];
1392: for (j=a->i[i]; j<a->i[i+1]; j++) {
1393: if (a->j[j] == i) {
1394: a->diag[i] = j;
1395: break;
1396: }
1397: }
1398: }
1399: return(0);
1400: }
1402: /*
1403: Checks for missing diagonals
1404: */
1407: PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d)
1408: {
1409: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1410: PetscInt *diag,*jj = a->j,i;
1413: *missing = PETSC_FALSE;
1414: if (A->rmap->n > 0 && !jj) {
1415: *missing = PETSC_TRUE;
1416: if (d) *d = 0;
1417: PetscInfo(A,"Matrix has no entries therefore is missing diagonal");
1418: } else {
1419: diag = a->diag;
1420: for (i=0; i<A->rmap->n; i++) {
1421: if (jj[diag[i]] != i) {
1422: *missing = PETSC_TRUE;
1423: if (d) *d = i;
1424: PetscInfo1(A,"Matrix is missing diagonal number %D",i);
1425: break;
1426: }
1427: }
1428: }
1429: return(0);
1430: }
1434: PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift)
1435: {
1436: Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data;
1438: PetscInt i,*diag,m = A->rmap->n;
1439: MatScalar *v = a->a;
1440: PetscScalar *idiag,*mdiag;
1443: if (a->idiagvalid) return(0);
1444: MatMarkDiagonal_SeqAIJ(A);
1445: diag = a->diag;
1446: if (!a->idiag) {
1447: PetscMalloc3(m,PetscScalar,&a->idiag,m,PetscScalar,&a->mdiag,m,PetscScalar,&a->ssor_work);
1448: PetscLogObjectMemory(A, 3*m*sizeof(PetscScalar));
1449: v = a->a;
1450: }
1451: mdiag = a->mdiag;
1452: idiag = a->idiag;
1454: if (omega == 1.0 && !PetscAbsScalar(fshift)) {
1455: for (i=0; i<m; i++) {
1456: mdiag[i] = v[diag[i]];
1457: if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i);
1458: idiag[i] = 1.0/v[diag[i]];
1459: }
1460: PetscLogFlops(m);
1461: } else {
1462: for (i=0; i<m; i++) {
1463: mdiag[i] = v[diag[i]];
1464: idiag[i] = omega/(fshift + v[diag[i]]);
1465: }
1466: PetscLogFlops(2.0*m);
1467: }
1468: a->idiagvalid = PETSC_TRUE;
1469: return(0);
1470: }
1472: #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h>
1475: PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
1476: {
1477: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1478: PetscScalar *x,d,sum,*t,scale;
1479: const MatScalar *v = a->a,*idiag=0,*mdiag;
1480: const PetscScalar *b, *bs,*xb, *ts;
1481: PetscErrorCode ierr;
1482: PetscInt n = A->cmap->n,m = A->rmap->n,i;
1483: const PetscInt *idx,*diag;
1486: its = its*lits;
1488: if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */
1489: if (!a->idiagvalid) {MatInvertDiagonal_SeqAIJ(A,omega,fshift);}
1490: a->fshift = fshift;
1491: a->omega = omega;
1493: diag = a->diag;
1494: t = a->ssor_work;
1495: idiag = a->idiag;
1496: mdiag = a->mdiag;
1498: VecGetArray(xx,&x);
1499: VecGetArrayRead(bb,&b);
1500: /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */
1501: if (flag == SOR_APPLY_UPPER) {
1502: /* apply (U + D/omega) to the vector */
1503: bs = b;
1504: for (i=0; i<m; i++) {
1505: d = fshift + mdiag[i];
1506: n = a->i[i+1] - diag[i] - 1;
1507: idx = a->j + diag[i] + 1;
1508: v = a->a + diag[i] + 1;
1509: sum = b[i]*d/omega;
1510: PetscSparseDensePlusDot(sum,bs,v,idx,n);
1511: x[i] = sum;
1512: }
1513: VecRestoreArray(xx,&x);
1514: VecRestoreArrayRead(bb,&b);
1515: PetscLogFlops(a->nz);
1516: return(0);
1517: }
1519: if (flag == SOR_APPLY_LOWER) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented");
1520: else if (flag & SOR_EISENSTAT) {
1521: /* Let A = L + U + D; where L is lower trianglar,
1522: U is upper triangular, E = D/omega; This routine applies
1524: (L + E)^{-1} A (U + E)^{-1}
1526: to a vector efficiently using Eisenstat's trick.
1527: */
1528: scale = (2.0/omega) - 1.0;
1530: /* x = (E + U)^{-1} b */
1531: for (i=m-1; i>=0; i--) {
1532: n = a->i[i+1] - diag[i] - 1;
1533: idx = a->j + diag[i] + 1;
1534: v = a->a + diag[i] + 1;
1535: sum = b[i];
1536: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1537: x[i] = sum*idiag[i];
1538: }
1540: /* t = b - (2*E - D)x */
1541: v = a->a;
1542: for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i];
1544: /* t = (E + L)^{-1}t */
1545: ts = t;
1546: diag = a->diag;
1547: for (i=0; i<m; i++) {
1548: n = diag[i] - a->i[i];
1549: idx = a->j + a->i[i];
1550: v = a->a + a->i[i];
1551: sum = t[i];
1552: PetscSparseDenseMinusDot(sum,ts,v,idx,n);
1553: t[i] = sum*idiag[i];
1554: /* x = x + t */
1555: x[i] += t[i];
1556: }
1558: PetscLogFlops(6.0*m-1 + 2.0*a->nz);
1559: VecRestoreArray(xx,&x);
1560: VecRestoreArrayRead(bb,&b);
1561: return(0);
1562: }
1563: if (flag & SOR_ZERO_INITIAL_GUESS) {
1564: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1565: for (i=0; i<m; i++) {
1566: n = diag[i] - a->i[i];
1567: idx = a->j + a->i[i];
1568: v = a->a + a->i[i];
1569: sum = b[i];
1570: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1571: t[i] = sum;
1572: x[i] = sum*idiag[i];
1573: }
1574: xb = t;
1575: PetscLogFlops(a->nz);
1576: } else xb = b;
1577: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1578: for (i=m-1; i>=0; i--) {
1579: n = a->i[i+1] - diag[i] - 1;
1580: idx = a->j + diag[i] + 1;
1581: v = a->a + diag[i] + 1;
1582: sum = xb[i];
1583: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1584: if (xb == b) {
1585: x[i] = sum*idiag[i];
1586: } else {
1587: x[i] = (1-omega)*x[i] + sum*idiag[i];
1588: }
1589: }
1590: PetscLogFlops(a->nz);
1591: }
1592: its--;
1593: }
1594: while (its--) {
1595: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1596: for (i=0; i<m; i++) {
1597: n = a->i[i+1] - a->i[i];
1598: idx = a->j + a->i[i];
1599: v = a->a + a->i[i];
1600: sum = b[i];
1601: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1602: x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i];
1603: }
1604: PetscLogFlops(2.0*a->nz);
1605: }
1606: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1607: for (i=m-1; i>=0; i--) {
1608: n = a->i[i+1] - a->i[i];
1609: idx = a->j + a->i[i];
1610: v = a->a + a->i[i];
1611: sum = b[i];
1612: PetscSparseDenseMinusDot(sum,x,v,idx,n);
1613: x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i];
1614: }
1615: PetscLogFlops(2.0*a->nz);
1616: }
1617: }
1618: VecRestoreArray(xx,&x);
1619: VecRestoreArrayRead(bb,&b);
1620: return(0);
1621: }
1626: PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info)
1627: {
1628: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1631: info->block_size = 1.0;
1632: info->nz_allocated = (double)a->maxnz;
1633: info->nz_used = (double)a->nz;
1634: info->nz_unneeded = (double)(a->maxnz - a->nz);
1635: info->assemblies = (double)A->num_ass;
1636: info->mallocs = (double)A->info.mallocs;
1637: info->memory = ((PetscObject)A)->mem;
1638: if (A->factortype) {
1639: info->fill_ratio_given = A->info.fill_ratio_given;
1640: info->fill_ratio_needed = A->info.fill_ratio_needed;
1641: info->factor_mallocs = A->info.factor_mallocs;
1642: } else {
1643: info->fill_ratio_given = 0;
1644: info->fill_ratio_needed = 0;
1645: info->factor_mallocs = 0;
1646: }
1647: return(0);
1648: }
1652: PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1653: {
1654: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1655: PetscInt i,m = A->rmap->n - 1,d = 0;
1656: PetscErrorCode ierr;
1657: const PetscScalar *xx;
1658: PetscScalar *bb;
1659: PetscBool missing;
1662: if (x && b) {
1663: VecGetArrayRead(x,&xx);
1664: VecGetArray(b,&bb);
1665: for (i=0; i<N; i++) {
1666: if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]);
1667: bb[rows[i]] = diag*xx[rows[i]];
1668: }
1669: VecRestoreArrayRead(x,&xx);
1670: VecRestoreArray(b,&bb);
1671: }
1673: if (a->keepnonzeropattern) {
1674: for (i=0; i<N; i++) {
1675: if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]);
1676: PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));
1677: }
1678: if (diag != 0.0) {
1679: MatMissingDiagonal_SeqAIJ(A,&missing,&d);
1680: if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d);
1681: for (i=0; i<N; i++) {
1682: a->a[a->diag[rows[i]]] = diag;
1683: }
1684: }
1685: A->same_nonzero = PETSC_TRUE;
1686: } else {
1687: if (diag != 0.0) {
1688: for (i=0; i<N; i++) {
1689: if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]);
1690: if (a->ilen[rows[i]] > 0) {
1691: a->ilen[rows[i]] = 1;
1692: a->a[a->i[rows[i]]] = diag;
1693: a->j[a->i[rows[i]]] = rows[i];
1694: } else { /* in case row was completely empty */
1695: MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);
1696: }
1697: }
1698: } else {
1699: for (i=0; i<N; i++) {
1700: if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]);
1701: a->ilen[rows[i]] = 0;
1702: }
1703: }
1704: A->same_nonzero = PETSC_FALSE;
1705: }
1706: MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);
1707: return(0);
1708: }
1712: PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1713: {
1714: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1715: PetscInt i,j,m = A->rmap->n - 1,d = 0;
1716: PetscErrorCode ierr;
1717: PetscBool missing,*zeroed,vecs = PETSC_FALSE;
1718: const PetscScalar *xx;
1719: PetscScalar *bb;
1722: if (x && b) {
1723: VecGetArrayRead(x,&xx);
1724: VecGetArray(b,&bb);
1725: vecs = PETSC_TRUE;
1726: }
1727: PetscMalloc(A->rmap->n*sizeof(PetscBool),&zeroed);
1728: PetscMemzero(zeroed,A->rmap->n*sizeof(PetscBool));
1729: for (i=0; i<N; i++) {
1730: if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]);
1731: PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));
1733: zeroed[rows[i]] = PETSC_TRUE;
1734: }
1735: for (i=0; i<A->rmap->n; i++) {
1736: if (!zeroed[i]) {
1737: for (j=a->i[i]; j<a->i[i+1]; j++) {
1738: if (zeroed[a->j[j]]) {
1739: if (vecs) bb[i] -= a->a[j]*xx[a->j[j]];
1740: a->a[j] = 0.0;
1741: }
1742: }
1743: } else if (vecs) bb[i] = diag*xx[i];
1744: }
1745: if (x && b) {
1746: VecRestoreArrayRead(x,&xx);
1747: VecRestoreArray(b,&bb);
1748: }
1749: PetscFree(zeroed);
1750: if (diag != 0.0) {
1751: MatMissingDiagonal_SeqAIJ(A,&missing,&d);
1752: if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d);
1753: for (i=0; i<N; i++) {
1754: a->a[a->diag[rows[i]]] = diag;
1755: }
1756: }
1757: A->same_nonzero = PETSC_TRUE;
1758: MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);
1759: return(0);
1760: }
1764: PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1765: {
1766: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1767: PetscInt *itmp;
1770: if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row);
1772: *nz = a->i[row+1] - a->i[row];
1773: if (v) *v = a->a + a->i[row];
1774: if (idx) {
1775: itmp = a->j + a->i[row];
1776: if (*nz) *idx = itmp;
1777: else *idx = 0;
1778: }
1779: return(0);
1780: }
1782: /* remove this function? */
1785: PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1786: {
1788: return(0);
1789: }
1793: PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm)
1794: {
1795: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1796: MatScalar *v = a->a;
1797: PetscReal sum = 0.0;
1799: PetscInt i,j;
1802: if (type == NORM_FROBENIUS) {
1803: for (i=0; i<a->nz; i++) {
1804: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1805: }
1806: *nrm = PetscSqrtReal(sum);
1807: } else if (type == NORM_1) {
1808: PetscReal *tmp;
1809: PetscInt *jj = a->j;
1810: PetscMalloc((A->cmap->n+1)*sizeof(PetscReal),&tmp);
1811: PetscMemzero(tmp,A->cmap->n*sizeof(PetscReal));
1812: *nrm = 0.0;
1813: for (j=0; j<a->nz; j++) {
1814: tmp[*jj++] += PetscAbsScalar(*v); v++;
1815: }
1816: for (j=0; j<A->cmap->n; j++) {
1817: if (tmp[j] > *nrm) *nrm = tmp[j];
1818: }
1819: PetscFree(tmp);
1820: } else if (type == NORM_INFINITY) {
1821: *nrm = 0.0;
1822: for (j=0; j<A->rmap->n; j++) {
1823: v = a->a + a->i[j];
1824: sum = 0.0;
1825: for (i=0; i<a->i[j+1]-a->i[j]; i++) {
1826: sum += PetscAbsScalar(*v); v++;
1827: }
1828: if (sum > *nrm) *nrm = sum;
1829: }
1830: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm");
1831: return(0);
1832: }
1834: /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */
1837: PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B)
1838: {
1840: PetscInt i,j,anzj;
1841: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b;
1842: PetscInt an=A->cmap->N,am=A->rmap->N;
1843: PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j;
1846: /* Allocate space for symbolic transpose info and work array */
1847: PetscMalloc((an+1)*sizeof(PetscInt),&ati);
1848: PetscMalloc(ai[am]*sizeof(PetscInt),&atj);
1849: PetscMalloc(an*sizeof(PetscInt),&atfill);
1850: PetscMemzero(ati,(an+1)*sizeof(PetscInt));
1852: /* Walk through aj and count ## of non-zeros in each row of A^T. */
1853: /* Note: offset by 1 for fast conversion into csr format. */
1854: for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1;
1855: /* Form ati for csr format of A^T. */
1856: for (i=0;i<an;i++) ati[i+1] += ati[i];
1858: /* Copy ati into atfill so we have locations of the next free space in atj */
1859: PetscMemcpy(atfill,ati,an*sizeof(PetscInt));
1861: /* Walk through A row-wise and mark nonzero entries of A^T. */
1862: for (i=0;i<am;i++) {
1863: anzj = ai[i+1] - ai[i];
1864: for (j=0;j<anzj;j++) {
1865: atj[atfill[*aj]] = i;
1866: atfill[*aj++] += 1;
1867: }
1868: }
1870: /* Clean up temporary space and complete requests. */
1871: PetscFree(atfill);
1872: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);
1874: (*B)->rmap->bs = A->cmap->bs;
1875: (*B)->cmap->bs = A->rmap->bs;
1877: b = (Mat_SeqAIJ*)((*B)->data);
1878: b->free_a = PETSC_FALSE;
1879: b->free_ij = PETSC_TRUE;
1880: b->nonew = 0;
1881: return(0);
1882: }
1886: PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B)
1887: {
1888: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1889: Mat C;
1891: PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col;
1892: MatScalar *array = a->a;
1895: if (reuse == MAT_REUSE_MATRIX && A == *B && m != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Square matrix only for in-place");
1897: if (reuse == MAT_INITIAL_MATRIX || *B == A) {
1898: PetscMalloc((1+A->cmap->n)*sizeof(PetscInt),&col);
1899: PetscMemzero(col,(1+A->cmap->n)*sizeof(PetscInt));
1901: for (i=0; i<ai[m]; i++) col[aj[i]] += 1;
1902: MatCreate(PetscObjectComm((PetscObject)A),&C);
1903: MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);
1904: MatSetBlockSizes(C,A->cmap->bs,A->rmap->bs);
1905: MatSetType(C,((PetscObject)A)->type_name);
1906: MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);
1907: PetscFree(col);
1908: } else {
1909: C = *B;
1910: }
1912: for (i=0; i<m; i++) {
1913: len = ai[i+1]-ai[i];
1914: MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);
1915: array += len;
1916: aj += len;
1917: }
1918: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1919: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1921: if (reuse == MAT_INITIAL_MATRIX || *B != A) {
1922: *B = C;
1923: } else {
1924: MatHeaderMerge(A,C);
1925: }
1926: return(0);
1927: }
1931: PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f)
1932: {
1933: Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) A->data;
1934: PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr;
1935: MatScalar *va,*vb;
1937: PetscInt ma,na,mb,nb, i;
1940: bij = (Mat_SeqAIJ*) B->data;
1942: MatGetSize(A,&ma,&na);
1943: MatGetSize(B,&mb,&nb);
1944: if (ma!=nb || na!=mb) {
1945: *f = PETSC_FALSE;
1946: return(0);
1947: }
1948: aii = aij->i; bii = bij->i;
1949: adx = aij->j; bdx = bij->j;
1950: va = aij->a; vb = bij->a;
1951: PetscMalloc(ma*sizeof(PetscInt),&aptr);
1952: PetscMalloc(mb*sizeof(PetscInt),&bptr);
1953: for (i=0; i<ma; i++) aptr[i] = aii[i];
1954: for (i=0; i<mb; i++) bptr[i] = bii[i];
1956: *f = PETSC_TRUE;
1957: for (i=0; i<ma; i++) {
1958: while (aptr[i]<aii[i+1]) {
1959: PetscInt idc,idr;
1960: PetscScalar vc,vr;
1961: /* column/row index/value */
1962: idc = adx[aptr[i]];
1963: idr = bdx[bptr[idc]];
1964: vc = va[aptr[i]];
1965: vr = vb[bptr[idc]];
1966: if (i!=idr || PetscAbsScalar(vc-vr) > tol) {
1967: *f = PETSC_FALSE;
1968: goto done;
1969: } else {
1970: aptr[i]++;
1971: if (B || i!=idc) bptr[idc]++;
1972: }
1973: }
1974: }
1975: done:
1976: PetscFree(aptr);
1977: PetscFree(bptr);
1978: return(0);
1979: }
1983: PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f)
1984: {
1985: Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) A->data;
1986: PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr;
1987: MatScalar *va,*vb;
1989: PetscInt ma,na,mb,nb, i;
1992: bij = (Mat_SeqAIJ*) B->data;
1994: MatGetSize(A,&ma,&na);
1995: MatGetSize(B,&mb,&nb);
1996: if (ma!=nb || na!=mb) {
1997: *f = PETSC_FALSE;
1998: return(0);
1999: }
2000: aii = aij->i; bii = bij->i;
2001: adx = aij->j; bdx = bij->j;
2002: va = aij->a; vb = bij->a;
2003: PetscMalloc(ma*sizeof(PetscInt),&aptr);
2004: PetscMalloc(mb*sizeof(PetscInt),&bptr);
2005: for (i=0; i<ma; i++) aptr[i] = aii[i];
2006: for (i=0; i<mb; i++) bptr[i] = bii[i];
2008: *f = PETSC_TRUE;
2009: for (i=0; i<ma; i++) {
2010: while (aptr[i]<aii[i+1]) {
2011: PetscInt idc,idr;
2012: PetscScalar vc,vr;
2013: /* column/row index/value */
2014: idc = adx[aptr[i]];
2015: idr = bdx[bptr[idc]];
2016: vc = va[aptr[i]];
2017: vr = vb[bptr[idc]];
2018: if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) {
2019: *f = PETSC_FALSE;
2020: goto done;
2021: } else {
2022: aptr[i]++;
2023: if (B || i!=idc) bptr[idc]++;
2024: }
2025: }
2026: }
2027: done:
2028: PetscFree(aptr);
2029: PetscFree(bptr);
2030: return(0);
2031: }
2035: PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f)
2036: {
2040: MatIsTranspose_SeqAIJ(A,A,tol,f);
2041: return(0);
2042: }
2046: PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f)
2047: {
2051: MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);
2052: return(0);
2053: }
2057: PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr)
2058: {
2059: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2060: PetscScalar *l,*r,x;
2061: MatScalar *v;
2063: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj;
2066: if (ll) {
2067: /* The local size is used so that VecMPI can be passed to this routine
2068: by MatDiagonalScale_MPIAIJ */
2069: VecGetLocalSize(ll,&m);
2070: if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length");
2071: VecGetArray(ll,&l);
2072: v = a->a;
2073: for (i=0; i<m; i++) {
2074: x = l[i];
2075: M = a->i[i+1] - a->i[i];
2076: for (j=0; j<M; j++) (*v++) *= x;
2077: }
2078: VecRestoreArray(ll,&l);
2079: PetscLogFlops(nz);
2080: }
2081: if (rr) {
2082: VecGetLocalSize(rr,&n);
2083: if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length");
2084: VecGetArray(rr,&r);
2085: v = a->a; jj = a->j;
2086: for (i=0; i<nz; i++) (*v++) *= r[*jj++];
2087: VecRestoreArray(rr,&r);
2088: PetscLogFlops(nz);
2089: }
2090: MatSeqAIJInvalidateDiagonal(A);
2091: return(0);
2092: }
2096: PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B)
2097: {
2098: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c;
2100: PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens;
2101: PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi;
2102: const PetscInt *irow,*icol;
2103: PetscInt nrows,ncols;
2104: PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen;
2105: MatScalar *a_new,*mat_a;
2106: Mat C;
2107: PetscBool stride,sorted;
2110: ISSorted(isrow,&sorted);
2111: if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted");
2112: ISSorted(iscol,&sorted);
2113: if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted");
2115: ISGetIndices(isrow,&irow);
2116: ISGetLocalSize(isrow,&nrows);
2117: ISGetLocalSize(iscol,&ncols);
2119: ISStrideGetInfo(iscol,&first,&step);
2120: PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);
2121: if (stride && step == 1) {
2122: /* special case of contiguous rows */
2123: PetscMalloc2(nrows,PetscInt,&lens,nrows,PetscInt,&starts);
2124: /* loop over new rows determining lens and starting points */
2125: for (i=0; i<nrows; i++) {
2126: kstart = ai[irow[i]];
2127: kend = kstart + ailen[irow[i]];
2128: for (k=kstart; k<kend; k++) {
2129: if (aj[k] >= first) {
2130: starts[i] = k;
2131: break;
2132: }
2133: }
2134: sum = 0;
2135: while (k < kend) {
2136: if (aj[k++] >= first+ncols) break;
2137: sum++;
2138: }
2139: lens[i] = sum;
2140: }
2141: /* create submatrix */
2142: if (scall == MAT_REUSE_MATRIX) {
2143: PetscInt n_cols,n_rows;
2144: MatGetSize(*B,&n_rows,&n_cols);
2145: if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size");
2146: MatZeroEntries(*B);
2147: C = *B;
2148: } else {
2149: PetscInt rbs,cbs;
2150: MatCreate(PetscObjectComm((PetscObject)A),&C);
2151: MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);
2152: ISGetBlockSize(isrow,&rbs);
2153: ISGetBlockSize(iscol,&cbs);
2154: MatSetBlockSizes(C,rbs,cbs);
2155: MatSetType(C,((PetscObject)A)->type_name);
2156: MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);
2157: }
2158: c = (Mat_SeqAIJ*)C->data;
2160: /* loop over rows inserting into submatrix */
2161: a_new = c->a;
2162: j_new = c->j;
2163: i_new = c->i;
2165: for (i=0; i<nrows; i++) {
2166: ii = starts[i];
2167: lensi = lens[i];
2168: for (k=0; k<lensi; k++) {
2169: *j_new++ = aj[ii+k] - first;
2170: }
2171: PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));
2172: a_new += lensi;
2173: i_new[i+1] = i_new[i] + lensi;
2174: c->ilen[i] = lensi;
2175: }
2176: PetscFree2(lens,starts);
2177: } else {
2178: ISGetIndices(iscol,&icol);
2179: PetscMalloc(oldcols*sizeof(PetscInt),&smap);
2180: PetscMemzero(smap,oldcols*sizeof(PetscInt));
2181: PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);
2182: for (i=0; i<ncols; i++) {
2183: #if defined(PETSC_USE_DEBUG)
2184: if (icol[i] >= oldcols) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requesting column beyond largest column icol[%D] %D <= A->cmap->n %D",i,icol[i],oldcols);
2185: #endif
2186: smap[icol[i]] = i+1;
2187: }
2189: /* determine lens of each row */
2190: for (i=0; i<nrows; i++) {
2191: kstart = ai[irow[i]];
2192: kend = kstart + a->ilen[irow[i]];
2193: lens[i] = 0;
2194: for (k=kstart; k<kend; k++) {
2195: if (smap[aj[k]]) {
2196: lens[i]++;
2197: }
2198: }
2199: }
2200: /* Create and fill new matrix */
2201: if (scall == MAT_REUSE_MATRIX) {
2202: PetscBool equal;
2204: c = (Mat_SeqAIJ*)((*B)->data);
2205: if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size");
2206: PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);
2207: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros");
2208: PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));
2209: C = *B;
2210: } else {
2211: PetscInt rbs,cbs;
2212: MatCreate(PetscObjectComm((PetscObject)A),&C);
2213: MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);
2214: ISGetBlockSize(isrow,&rbs);
2215: ISGetBlockSize(iscol,&cbs);
2216: MatSetBlockSizes(C,rbs,cbs);
2217: MatSetType(C,((PetscObject)A)->type_name);
2218: MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);
2219: }
2220: c = (Mat_SeqAIJ*)(C->data);
2221: for (i=0; i<nrows; i++) {
2222: row = irow[i];
2223: kstart = ai[row];
2224: kend = kstart + a->ilen[row];
2225: mat_i = c->i[i];
2226: mat_j = c->j + mat_i;
2227: mat_a = c->a + mat_i;
2228: mat_ilen = c->ilen + i;
2229: for (k=kstart; k<kend; k++) {
2230: if ((tcol=smap[a->j[k]])) {
2231: *mat_j++ = tcol - 1;
2232: *mat_a++ = a->a[k];
2233: (*mat_ilen)++;
2235: }
2236: }
2237: }
2238: /* Free work space */
2239: ISRestoreIndices(iscol,&icol);
2240: PetscFree(smap);
2241: PetscFree(lens);
2242: }
2243: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2244: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2246: ISRestoreIndices(isrow,&irow);
2247: *B = C;
2248: return(0);
2249: }
2253: PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat)
2254: {
2256: Mat B;
2259: MatCreate(subComm,&B);
2260: MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);
2261: MatSetBlockSizes(B,mat->rmap->bs,mat->cmap->bs);
2262: MatSetType(B,MATSEQAIJ);
2263: MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);
2264: *subMat = B;
2265: return(0);
2266: }
2270: PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info)
2271: {
2272: Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data;
2274: Mat outA;
2275: PetscBool row_identity,col_identity;
2278: if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu");
2280: ISIdentity(row,&row_identity);
2281: ISIdentity(col,&col_identity);
2283: outA = inA;
2284: outA->factortype = MAT_FACTOR_LU;
2286: PetscObjectReference((PetscObject)row);
2287: ISDestroy(&a->row);
2289: a->row = row;
2291: PetscObjectReference((PetscObject)col);
2292: ISDestroy(&a->col);
2294: a->col = col;
2296: /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */
2297: ISDestroy(&a->icol);
2298: ISInvertPermutation(col,PETSC_DECIDE,&a->icol);
2299: PetscLogObjectParent(inA,a->icol);
2301: if (!a->solve_work) { /* this matrix may have been factored before */
2302: PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);
2303: PetscLogObjectMemory(inA, (inA->rmap->n+1)*sizeof(PetscScalar));
2304: }
2306: MatMarkDiagonal_SeqAIJ(inA);
2307: if (row_identity && col_identity) {
2308: MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);
2309: } else {
2310: MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);
2311: }
2312: return(0);
2313: }
2317: PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha)
2318: {
2319: Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data;
2320: PetscScalar oalpha = alpha;
2322: PetscBLASInt one = 1,bnz;
2325: PetscBLASIntCast(a->nz,&bnz);
2326: PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one));
2327: PetscLogFlops(a->nz);
2328: MatSeqAIJInvalidateDiagonal(inA);
2329: return(0);
2330: }
2334: PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
2335: {
2337: PetscInt i;
2340: if (scall == MAT_INITIAL_MATRIX) {
2341: PetscMalloc((n+1)*sizeof(Mat),B);
2342: }
2344: for (i=0; i<n; i++) {
2345: MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
2346: }
2347: return(0);
2348: }
2352: PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov)
2353: {
2354: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2356: PetscInt row,i,j,k,l,m,n,*nidx,isz,val;
2357: const PetscInt *idx;
2358: PetscInt start,end,*ai,*aj;
2359: PetscBT table;
2362: m = A->rmap->n;
2363: ai = a->i;
2364: aj = a->j;
2366: if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used");
2368: PetscMalloc((m+1)*sizeof(PetscInt),&nidx);
2369: PetscBTCreate(m,&table);
2371: for (i=0; i<is_max; i++) {
2372: /* Initialize the two local arrays */
2373: isz = 0;
2374: PetscBTMemzero(m,table);
2376: /* Extract the indices, assume there can be duplicate entries */
2377: ISGetIndices(is[i],&idx);
2378: ISGetLocalSize(is[i],&n);
2380: /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */
2381: for (j=0; j<n; ++j) {
2382: if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j];
2383: }
2384: ISRestoreIndices(is[i],&idx);
2385: ISDestroy(&is[i]);
2387: k = 0;
2388: for (j=0; j<ov; j++) { /* for each overlap */
2389: n = isz;
2390: for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */
2391: row = nidx[k];
2392: start = ai[row];
2393: end = ai[row+1];
2394: for (l = start; l<end; l++) {
2395: val = aj[l];
2396: if (!PetscBTLookupSet(table,val)) nidx[isz++] = val;
2397: }
2398: }
2399: }
2400: ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));
2401: }
2402: PetscBTDestroy(&table);
2403: PetscFree(nidx);
2404: return(0);
2405: }
2407: /* -------------------------------------------------------------- */
2410: PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B)
2411: {
2412: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2414: PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n;
2415: const PetscInt *row,*col;
2416: PetscInt *cnew,j,*lens;
2417: IS icolp,irowp;
2418: PetscInt *cwork = NULL;
2419: PetscScalar *vwork = NULL;
2422: ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);
2423: ISGetIndices(irowp,&row);
2424: ISInvertPermutation(colp,PETSC_DECIDE,&icolp);
2425: ISGetIndices(icolp,&col);
2427: /* determine lengths of permuted rows */
2428: PetscMalloc((m+1)*sizeof(PetscInt),&lens);
2429: for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i];
2430: MatCreate(PetscObjectComm((PetscObject)A),B);
2431: MatSetSizes(*B,m,n,m,n);
2432: MatSetBlockSizes(*B,A->rmap->bs,A->cmap->bs);
2433: MatSetType(*B,((PetscObject)A)->type_name);
2434: MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);
2435: PetscFree(lens);
2437: PetscMalloc(n*sizeof(PetscInt),&cnew);
2438: for (i=0; i<m; i++) {
2439: MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);
2440: for (j=0; j<nz; j++) cnew[j] = col[cwork[j]];
2441: MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);
2442: MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);
2443: }
2444: PetscFree(cnew);
2446: (*B)->assembled = PETSC_FALSE;
2448: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
2449: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
2450: ISRestoreIndices(irowp,&row);
2451: ISRestoreIndices(icolp,&col);
2452: ISDestroy(&irowp);
2453: ISDestroy(&icolp);
2454: return(0);
2455: }
2459: PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str)
2460: {
2464: /* If the two matrices have the same copy implementation, use fast copy. */
2465: if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) {
2466: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2467: Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;
2469: if (a->i[A->rmap->n] != b->i[B->rmap->n]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different");
2470: PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));
2471: } else {
2472: MatCopy_Basic(A,B,str);
2473: }
2474: return(0);
2475: }
2479: PetscErrorCode MatSetUp_SeqAIJ(Mat A)
2480: {
2484: MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);
2485: return(0);
2486: }
2490: PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[])
2491: {
2492: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2495: *array = a->a;
2496: return(0);
2497: }
2501: PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[])
2502: {
2504: return(0);
2505: }
2509: PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx)
2510: {
2511: PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void*))coloring->f;
2513: PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow;
2514: PetscScalar dx,*y,*xx,*w3_array;
2515: PetscScalar *vscale_array;
2516: PetscReal epsilon = coloring->error_rel,umin = coloring->umin;
2517: Vec w1,w2,w3;
2518: void *fctx = coloring->fctx;
2519: PetscBool flg = PETSC_FALSE;
2522: if (!coloring->w1) {
2523: VecDuplicate(x1,&coloring->w1);
2524: PetscLogObjectParent(coloring,coloring->w1);
2525: VecDuplicate(x1,&coloring->w2);
2526: PetscLogObjectParent(coloring,coloring->w2);
2527: VecDuplicate(x1,&coloring->w3);
2528: PetscLogObjectParent(coloring,coloring->w3);
2529: }
2530: w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3;
2532: MatSetUnfactored(J);
2533: PetscOptionsGetBool(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,NULL);
2534: if (flg) {
2535: PetscInfo(coloring,"Not calling MatZeroEntries()\n");
2536: } else {
2537: PetscBool assembled;
2538: MatAssembled(J,&assembled);
2539: if (assembled) {
2540: MatZeroEntries(J);
2541: }
2542: }
2544: VecGetOwnershipRange(x1,&start,&end);
2545: VecGetSize(x1,&N);
2547: if (!coloring->fset) {
2548: PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
2549: (*f)(sctx,x1,w1,fctx);
2550: PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
2551: } else {
2552: coloring->fset = PETSC_FALSE;
2553: }
2555: /*
2556: Compute all the scale factors and share with other processors
2557: */
2558: VecGetArray(x1,&xx);xx = xx - start;
2559: VecGetArray(coloring->vscale,&vscale_array);vscale_array = vscale_array - start;
2560: for (k=0; k<coloring->ncolors; k++) {
2561: /*
2562: Loop over each column associated with color adding the
2563: perturbation to the vector w3.
2564: */
2565: for (l=0; l<coloring->ncolumns[k]; l++) {
2566: col = coloring->columns[k][l]; /* column of the matrix we are probing for */
2567: dx = xx[col];
2568: if (dx == 0.0) dx = 1.0;
2569: if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin;
2570: else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
2571: dx *= epsilon;
2572: vscale_array[col] = 1.0/dx;
2573: }
2574: }
2575: vscale_array = vscale_array + start;
2577: VecRestoreArray(coloring->vscale,&vscale_array);
2578: VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
2579: VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
2581: /* VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD);
2582: VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/
2584: if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow;
2585: else vscaleforrow = coloring->columnsforrow;
2587: VecGetArray(coloring->vscale,&vscale_array);
2588: /*
2589: Loop over each color
2590: */
2591: for (k=0; k<coloring->ncolors; k++) {
2592: coloring->currentcolor = k;
2594: VecCopy(x1,w3);
2595: VecGetArray(w3,&w3_array);w3_array = w3_array - start;
2596: /*
2597: Loop over each column associated with color adding the
2598: perturbation to the vector w3.
2599: */
2600: for (l=0; l<coloring->ncolumns[k]; l++) {
2601: col = coloring->columns[k][l]; /* column of the matrix we are probing for */
2602: dx = xx[col];
2603: if (dx == 0.0) dx = 1.0;
2604: if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin;
2605: else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
2606: dx *= epsilon;
2607: if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter");
2608: w3_array[col] += dx;
2609: }
2610: w3_array = w3_array + start; VecRestoreArray(w3,&w3_array);
2612: /*
2613: Evaluate function at x1 + dx (here dx is a vector of perturbations)
2614: */
2616: PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
2617: (*f)(sctx,w3,w2,fctx);
2618: PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
2619: VecAXPY(w2,-1.0,w1);
2621: /*
2622: Loop over rows of vector, putting results into Jacobian matrix
2623: */
2624: VecGetArray(w2,&y);
2625: for (l=0; l<coloring->nrows[k]; l++) {
2626: row = coloring->rows[k][l];
2627: col = coloring->columnsforrow[k][l];
2628: y[row] *= vscale_array[vscaleforrow[k][l]];
2629: srow = row + start;
2630: MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);
2631: }
2632: VecRestoreArray(w2,&y);
2633: }
2634: coloring->currentcolor = k;
2636: VecRestoreArray(coloring->vscale,&vscale_array);
2637: xx = xx + start;
2638: VecRestoreArray(x1,&xx);
2639: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
2640: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
2641: return(0);
2642: }
2644: /*
2645: Computes the number of nonzeros per row needed for preallocation when X and Y
2646: have different nonzero structure.
2647: */
2650: PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz)
2651: {
2652: PetscInt i,m=Y->rmap->N;
2653: Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data;
2654: Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data;
2655: const PetscInt *xi = x->i,*yi = y->i;
2658: /* Set the number of nonzeros in the new matrix */
2659: for (i=0; i<m; i++) {
2660: PetscInt j,k,nzx = xi[i+1] - xi[i],nzy = yi[i+1] - yi[i];
2661: const PetscInt *xj = x->j+xi[i],*yj = y->j+yi[i];
2662: nnz[i] = 0;
2663: for (j=0,k=0; j<nzx; j++) { /* Point in X */
2664: for (; k<nzy && yj[k]<xj[j]; k++) nnz[i]++; /* Catch up to X */
2665: if (k<nzy && yj[k]==xj[j]) k++; /* Skip duplicate */
2666: nnz[i]++;
2667: }
2668: for (; k<nzy; k++) nnz[i]++;
2669: }
2670: return(0);
2671: }
2675: PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)
2676: {
2678: PetscInt i;
2679: Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data;
2680: PetscBLASInt one=1,bnz;
2683: PetscBLASIntCast(x->nz,&bnz);
2684: if (str == SAME_NONZERO_PATTERN) {
2685: PetscScalar alpha = a;
2686: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one));
2687: MatSeqAIJInvalidateDiagonal(Y);
2688: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
2689: if (y->xtoy && y->XtoY != X) {
2690: PetscFree(y->xtoy);
2691: MatDestroy(&y->XtoY);
2692: }
2693: if (!y->xtoy) { /* get xtoy */
2694: MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,NULL, y->i,y->j,NULL, &y->xtoy);
2695: y->XtoY = X;
2696: PetscObjectReference((PetscObject)X);
2697: }
2698: for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]);
2699: PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %g\n",x->nz,y->nz,(double)(PetscReal)(x->nz)/(y->nz+1));
2700: } else {
2701: Mat B;
2702: PetscInt *nnz;
2703: PetscMalloc(Y->rmap->N*sizeof(PetscInt),&nnz);
2704: MatCreate(PetscObjectComm((PetscObject)Y),&B);
2705: PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);
2706: MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);
2707: MatSetBlockSizes(B,Y->rmap->bs,Y->cmap->bs);
2708: MatSetType(B,(MatType) ((PetscObject)Y)->type_name);
2709: MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);
2710: MatSeqAIJSetPreallocation(B,0,nnz);
2711: MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
2712: MatHeaderReplace(Y,B);
2713: PetscFree(nnz);
2714: }
2715: return(0);
2716: }
2720: PetscErrorCode MatConjugate_SeqAIJ(Mat mat)
2721: {
2722: #if defined(PETSC_USE_COMPLEX)
2723: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
2724: PetscInt i,nz;
2725: PetscScalar *a;
2728: nz = aij->nz;
2729: a = aij->a;
2730: for (i=0; i<nz; i++) a[i] = PetscConj(a[i]);
2731: #else
2733: #endif
2734: return(0);
2735: }
2739: PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[])
2740: {
2741: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2743: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
2744: PetscReal atmp;
2745: PetscScalar *x;
2746: MatScalar *aa;
2749: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2750: aa = a->a;
2751: ai = a->i;
2752: aj = a->j;
2754: VecSet(v,0.0);
2755: VecGetArray(v,&x);
2756: VecGetLocalSize(v,&n);
2757: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2758: for (i=0; i<m; i++) {
2759: ncols = ai[1] - ai[0]; ai++;
2760: x[i] = 0.0;
2761: for (j=0; j<ncols; j++) {
2762: atmp = PetscAbsScalar(*aa);
2763: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;}
2764: aa++; aj++;
2765: }
2766: }
2767: VecRestoreArray(v,&x);
2768: return(0);
2769: }
2773: PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[])
2774: {
2775: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2777: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
2778: PetscScalar *x;
2779: MatScalar *aa;
2782: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2783: aa = a->a;
2784: ai = a->i;
2785: aj = a->j;
2787: VecSet(v,0.0);
2788: VecGetArray(v,&x);
2789: VecGetLocalSize(v,&n);
2790: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2791: for (i=0; i<m; i++) {
2792: ncols = ai[1] - ai[0]; ai++;
2793: if (ncols == A->cmap->n) { /* row is dense */
2794: x[i] = *aa; if (idx) idx[i] = 0;
2795: } else { /* row is sparse so already KNOW maximum is 0.0 or higher */
2796: x[i] = 0.0;
2797: if (idx) {
2798: idx[i] = 0; /* in case ncols is zero */
2799: for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */
2800: if (aj[j] > j) {
2801: idx[i] = j;
2802: break;
2803: }
2804: }
2805: }
2806: }
2807: for (j=0; j<ncols; j++) {
2808: if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;}
2809: aa++; aj++;
2810: }
2811: }
2812: VecRestoreArray(v,&x);
2813: return(0);
2814: }
2818: PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[])
2819: {
2820: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2822: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
2823: PetscReal atmp;
2824: PetscScalar *x;
2825: MatScalar *aa;
2828: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2829: aa = a->a;
2830: ai = a->i;
2831: aj = a->j;
2833: VecSet(v,0.0);
2834: VecGetArray(v,&x);
2835: VecGetLocalSize(v,&n);
2836: if (n != A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector, %d vs. %d rows", A->rmap->n, n);
2837: for (i=0; i<m; i++) {
2838: ncols = ai[1] - ai[0]; ai++;
2839: if (ncols) {
2840: /* Get first nonzero */
2841: for (j = 0; j < ncols; j++) {
2842: atmp = PetscAbsScalar(aa[j]);
2843: if (atmp > 1.0e-12) {
2844: x[i] = atmp;
2845: if (idx) idx[i] = aj[j];
2846: break;
2847: }
2848: }
2849: if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;}
2850: } else {
2851: x[i] = 0.0; if (idx) idx[i] = 0;
2852: }
2853: for (j = 0; j < ncols; j++) {
2854: atmp = PetscAbsScalar(*aa);
2855: if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;}
2856: aa++; aj++;
2857: }
2858: }
2859: VecRestoreArray(v,&x);
2860: return(0);
2861: }
2865: PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[])
2866: {
2867: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2869: PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n;
2870: PetscScalar *x;
2871: MatScalar *aa;
2874: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2875: aa = a->a;
2876: ai = a->i;
2877: aj = a->j;
2879: VecSet(v,0.0);
2880: VecGetArray(v,&x);
2881: VecGetLocalSize(v,&n);
2882: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2883: for (i=0; i<m; i++) {
2884: ncols = ai[1] - ai[0]; ai++;
2885: if (ncols == A->cmap->n) { /* row is dense */
2886: x[i] = *aa; if (idx) idx[i] = 0;
2887: } else { /* row is sparse so already KNOW minimum is 0.0 or lower */
2888: x[i] = 0.0;
2889: if (idx) { /* find first implicit 0.0 in the row */
2890: idx[i] = 0; /* in case ncols is zero */
2891: for (j=0; j<ncols; j++) {
2892: if (aj[j] > j) {
2893: idx[i] = j;
2894: break;
2895: }
2896: }
2897: }
2898: }
2899: for (j=0; j<ncols; j++) {
2900: if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;}
2901: aa++; aj++;
2902: }
2903: }
2904: VecRestoreArray(v,&x);
2905: return(0);
2906: }
2908: #include <petscblaslapack.h>
2909: #include <petsc-private/kernels/blockinvert.h>
2913: PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values)
2914: {
2915: Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data;
2917: PetscInt i,bs = A->rmap->bs,mbs = A->rmap->n/A->rmap->bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j;
2918: MatScalar *diag,work[25],*v_work;
2919: PetscReal shift = 0.0;
2922: if (a->ibdiagvalid) {
2923: if (values) *values = a->ibdiag;
2924: return(0);
2925: }
2926: MatMarkDiagonal_SeqAIJ(A);
2927: if (!a->ibdiag) {
2928: PetscMalloc(bs2*mbs*sizeof(PetscScalar),&a->ibdiag);
2929: PetscLogObjectMemory(A,bs2*mbs*sizeof(PetscScalar));
2930: }
2931: diag = a->ibdiag;
2932: if (values) *values = a->ibdiag;
2933: /* factor and invert each block */
2934: switch (bs) {
2935: case 1:
2936: for (i=0; i<mbs; i++) {
2937: MatGetValues(A,1,&i,1,&i,diag+i);
2938: diag[i] = (PetscScalar)1.0 / (diag[i] + shift);
2939: }
2940: break;
2941: case 2:
2942: for (i=0; i<mbs; i++) {
2943: ij[0] = 2*i; ij[1] = 2*i + 1;
2944: MatGetValues(A,2,ij,2,ij,diag);
2945: PetscKernel_A_gets_inverse_A_2(diag,shift);
2946: PetscKernel_A_gets_transpose_A_2(diag);
2947: diag += 4;
2948: }
2949: break;
2950: case 3:
2951: for (i=0; i<mbs; i++) {
2952: ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2;
2953: MatGetValues(A,3,ij,3,ij,diag);
2954: PetscKernel_A_gets_inverse_A_3(diag,shift);
2955: PetscKernel_A_gets_transpose_A_3(diag);
2956: diag += 9;
2957: }
2958: break;
2959: case 4:
2960: for (i=0; i<mbs; i++) {
2961: ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3;
2962: MatGetValues(A,4,ij,4,ij,diag);
2963: PetscKernel_A_gets_inverse_A_4(diag,shift);
2964: PetscKernel_A_gets_transpose_A_4(diag);
2965: diag += 16;
2966: }
2967: break;
2968: case 5:
2969: for (i=0; i<mbs; i++) {
2970: ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4;
2971: MatGetValues(A,5,ij,5,ij,diag);
2972: PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift);
2973: PetscKernel_A_gets_transpose_A_5(diag);
2974: diag += 25;
2975: }
2976: break;
2977: case 6:
2978: for (i=0; i<mbs; i++) {
2979: ij[0] = 6*i; ij[1] = 6*i + 1; ij[2] = 6*i + 2; ij[3] = 6*i + 3; ij[4] = 6*i + 4; ij[5] = 6*i + 5;
2980: MatGetValues(A,6,ij,6,ij,diag);
2981: PetscKernel_A_gets_inverse_A_6(diag,shift);
2982: PetscKernel_A_gets_transpose_A_6(diag);
2983: diag += 36;
2984: }
2985: break;
2986: case 7:
2987: for (i=0; i<mbs; i++) {
2988: ij[0] = 7*i; ij[1] = 7*i + 1; ij[2] = 7*i + 2; ij[3] = 7*i + 3; ij[4] = 7*i + 4; ij[5] = 7*i + 5; ij[5] = 7*i + 6;
2989: MatGetValues(A,7,ij,7,ij,diag);
2990: PetscKernel_A_gets_inverse_A_7(diag,shift);
2991: PetscKernel_A_gets_transpose_A_7(diag);
2992: diag += 49;
2993: }
2994: break;
2995: default:
2996: PetscMalloc3(bs,MatScalar,&v_work,bs,PetscInt,&v_pivots,bs,PetscInt,&IJ);
2997: for (i=0; i<mbs; i++) {
2998: for (j=0; j<bs; j++) {
2999: IJ[j] = bs*i + j;
3000: }
3001: MatGetValues(A,bs,IJ,bs,IJ,diag);
3002: PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work);
3003: PetscKernel_A_gets_transpose_A_N(diag,bs);
3004: diag += bs2;
3005: }
3006: PetscFree3(v_work,v_pivots,IJ);
3007: }
3008: a->ibdiagvalid = PETSC_TRUE;
3009: return(0);
3010: }
3014: static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx)
3015: {
3017: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data;
3018: PetscScalar a;
3019: PetscInt m,n,i,j,col;
3022: if (!x->assembled) {
3023: MatGetSize(x,&m,&n);
3024: for (i=0; i<m; i++) {
3025: for (j=0; j<aij->imax[i]; j++) {
3026: PetscRandomGetValue(rctx,&a);
3027: col = (PetscInt)(n*PetscRealPart(a));
3028: MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);
3029: }
3030: }
3031: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded");
3032: MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);
3033: MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);
3034: return(0);
3035: }
3037: extern PetscErrorCode MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*);
3038: /* -------------------------------------------------------------------*/
3039: static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ,
3040: MatGetRow_SeqAIJ,
3041: MatRestoreRow_SeqAIJ,
3042: MatMult_SeqAIJ,
3043: /* 4*/ MatMultAdd_SeqAIJ,
3044: MatMultTranspose_SeqAIJ,
3045: MatMultTransposeAdd_SeqAIJ,
3046: 0,
3047: 0,
3048: 0,
3049: /* 10*/ 0,
3050: MatLUFactor_SeqAIJ,
3051: 0,
3052: MatSOR_SeqAIJ,
3053: MatTranspose_SeqAIJ,
3054: /*1 5*/ MatGetInfo_SeqAIJ,
3055: MatEqual_SeqAIJ,
3056: MatGetDiagonal_SeqAIJ,
3057: MatDiagonalScale_SeqAIJ,
3058: MatNorm_SeqAIJ,
3059: /* 20*/ 0,
3060: MatAssemblyEnd_SeqAIJ,
3061: MatSetOption_SeqAIJ,
3062: MatZeroEntries_SeqAIJ,
3063: /* 24*/ MatZeroRows_SeqAIJ,
3064: 0,
3065: 0,
3066: 0,
3067: 0,
3068: /* 29*/ MatSetUp_SeqAIJ,
3069: 0,
3070: 0,
3071: 0,
3072: 0,
3073: /* 34*/ MatDuplicate_SeqAIJ,
3074: 0,
3075: 0,
3076: MatILUFactor_SeqAIJ,
3077: 0,
3078: /* 39*/ MatAXPY_SeqAIJ,
3079: MatGetSubMatrices_SeqAIJ,
3080: MatIncreaseOverlap_SeqAIJ,
3081: MatGetValues_SeqAIJ,
3082: MatCopy_SeqAIJ,
3083: /* 44*/ MatGetRowMax_SeqAIJ,
3084: MatScale_SeqAIJ,
3085: 0,
3086: MatDiagonalSet_SeqAIJ,
3087: MatZeroRowsColumns_SeqAIJ,
3088: /* 49*/ MatSetRandom_SeqAIJ,
3089: MatGetRowIJ_SeqAIJ,
3090: MatRestoreRowIJ_SeqAIJ,
3091: MatGetColumnIJ_SeqAIJ,
3092: MatRestoreColumnIJ_SeqAIJ,
3093: /* 54*/ MatFDColoringCreate_SeqAIJ,
3094: 0,
3095: 0,
3096: MatPermute_SeqAIJ,
3097: 0,
3098: /* 59*/ 0,
3099: MatDestroy_SeqAIJ,
3100: MatView_SeqAIJ,
3101: 0,
3102: MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ,
3103: /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ,
3104: MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ,
3105: 0,
3106: 0,
3107: 0,
3108: /* 69*/ MatGetRowMaxAbs_SeqAIJ,
3109: MatGetRowMinAbs_SeqAIJ,
3110: 0,
3111: MatSetColoring_SeqAIJ,
3112: 0,
3113: /* 74*/ MatSetValuesAdifor_SeqAIJ,
3114: MatFDColoringApply_AIJ,
3115: 0,
3116: 0,
3117: 0,
3118: /* 79*/ MatFindZeroDiagonals_SeqAIJ,
3119: 0,
3120: 0,
3121: 0,
3122: MatLoad_SeqAIJ,
3123: /* 84*/ MatIsSymmetric_SeqAIJ,
3124: MatIsHermitian_SeqAIJ,
3125: 0,
3126: 0,
3127: 0,
3128: /* 89*/ MatMatMult_SeqAIJ_SeqAIJ,
3129: MatMatMultSymbolic_SeqAIJ_SeqAIJ,
3130: MatMatMultNumeric_SeqAIJ_SeqAIJ,
3131: MatPtAP_SeqAIJ_SeqAIJ,
3132: MatPtAPSymbolic_SeqAIJ_SeqAIJ,
3133: /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ,
3134: MatMatTransposeMult_SeqAIJ_SeqAIJ,
3135: MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ,
3136: MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ,
3137: 0,
3138: /* 99*/ 0,
3139: 0,
3140: 0,
3141: MatConjugate_SeqAIJ,
3142: 0,
3143: /*104*/ MatSetValuesRow_SeqAIJ,
3144: MatRealPart_SeqAIJ,
3145: MatImaginaryPart_SeqAIJ,
3146: 0,
3147: 0,
3148: /*109*/ MatMatSolve_SeqAIJ,
3149: 0,
3150: MatGetRowMin_SeqAIJ,
3151: 0,
3152: MatMissingDiagonal_SeqAIJ,
3153: /*114*/ 0,
3154: 0,
3155: 0,
3156: 0,
3157: 0,
3158: /*119*/ 0,
3159: 0,
3160: 0,
3161: 0,
3162: MatGetMultiProcBlock_SeqAIJ,
3163: /*124*/ MatFindNonzeroRows_SeqAIJ,
3164: MatGetColumnNorms_SeqAIJ,
3165: MatInvertBlockDiagonal_SeqAIJ,
3166: 0,
3167: 0,
3168: /*129*/ 0,
3169: MatTransposeMatMult_SeqAIJ_SeqAIJ,
3170: MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ,
3171: MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ,
3172: MatTransposeColoringCreate_SeqAIJ,
3173: /*134*/ MatTransColoringApplySpToDen_SeqAIJ,
3174: MatTransColoringApplyDenToSp_SeqAIJ,
3175: MatRARt_SeqAIJ_SeqAIJ,
3176: MatRARtSymbolic_SeqAIJ_SeqAIJ,
3177: MatRARtNumeric_SeqAIJ_SeqAIJ,
3178: /*139*/0,
3179: 0
3180: };
3184: PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices)
3185: {
3186: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3187: PetscInt i,nz,n;
3190: nz = aij->maxnz;
3191: n = mat->rmap->n;
3192: for (i=0; i<nz; i++) {
3193: aij->j[i] = indices[i];
3194: }
3195: aij->nz = nz;
3196: for (i=0; i<n; i++) {
3197: aij->ilen[i] = aij->imax[i];
3198: }
3199: return(0);
3200: }
3204: /*@
3205: MatSeqAIJSetColumnIndices - Set the column indices for all the rows
3206: in the matrix.
3208: Input Parameters:
3209: + mat - the SeqAIJ matrix
3210: - indices - the column indices
3212: Level: advanced
3214: Notes:
3215: This can be called if you have precomputed the nonzero structure of the
3216: matrix and want to provide it to the matrix object to improve the performance
3217: of the MatSetValues() operation.
3219: You MUST have set the correct numbers of nonzeros per row in the call to
3220: MatCreateSeqAIJ(), and the columns indices MUST be sorted.
3222: MUST be called before any calls to MatSetValues();
3224: The indices should start with zero, not one.
3226: @*/
3227: PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices)
3228: {
3234: PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));
3235: return(0);
3236: }
3238: /* ----------------------------------------------------------------------------------------*/
3242: PetscErrorCode MatStoreValues_SeqAIJ(Mat mat)
3243: {
3244: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3246: size_t nz = aij->i[mat->rmap->n];
3249: if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
3251: /* allocate space for values if not already there */
3252: if (!aij->saved_values) {
3253: PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);
3254: PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));
3255: }
3257: /* copy values over */
3258: PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));
3259: return(0);
3260: }
3264: /*@
3265: MatStoreValues - Stashes a copy of the matrix values; this allows, for
3266: example, reuse of the linear part of a Jacobian, while recomputing the
3267: nonlinear portion.
3269: Collect on Mat
3271: Input Parameters:
3272: . mat - the matrix (currently only AIJ matrices support this option)
3274: Level: advanced
3276: Common Usage, with SNESSolve():
3277: $ Create Jacobian matrix
3278: $ Set linear terms into matrix
3279: $ Apply boundary conditions to matrix, at this time matrix must have
3280: $ final nonzero structure (i.e. setting the nonlinear terms and applying
3281: $ boundary conditions again will not change the nonzero structure
3282: $ MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);
3283: $ MatStoreValues(mat);
3284: $ Call SNESSetJacobian() with matrix
3285: $ In your Jacobian routine
3286: $ MatRetrieveValues(mat);
3287: $ Set nonlinear terms in matrix
3289: Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself:
3290: $ // build linear portion of Jacobian
3291: $ MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);
3292: $ MatStoreValues(mat);
3293: $ loop over nonlinear iterations
3294: $ MatRetrieveValues(mat);
3295: $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian
3296: $ // call MatAssemblyBegin/End() on matrix
3297: $ Solve linear system with Jacobian
3298: $ endloop
3300: Notes:
3301: Matrix must already be assemblied before calling this routine
3302: Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before
3303: calling this routine.
3305: When this is called multiple times it overwrites the previous set of stored values
3306: and does not allocated additional space.
3308: .seealso: MatRetrieveValues()
3310: @*/
3311: PetscErrorCode MatStoreValues(Mat mat)
3312: {
3317: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3318: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3319: PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));
3320: return(0);
3321: }
3325: PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat)
3326: {
3327: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
3329: PetscInt nz = aij->i[mat->rmap->n];
3332: if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
3333: if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first");
3334: /* copy values over */
3335: PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));
3336: return(0);
3337: }
3341: /*@
3342: MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for
3343: example, reuse of the linear part of a Jacobian, while recomputing the
3344: nonlinear portion.
3346: Collect on Mat
3348: Input Parameters:
3349: . mat - the matrix (currently on AIJ matrices support this option)
3351: Level: advanced
3353: .seealso: MatStoreValues()
3355: @*/
3356: PetscErrorCode MatRetrieveValues(Mat mat)
3357: {
3362: if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3363: if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3364: PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));
3365: return(0);
3366: }
3369: /* --------------------------------------------------------------------------------*/
3372: /*@C
3373: MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format
3374: (the default parallel PETSc format). For good matrix assembly performance
3375: the user should preallocate the matrix storage by setting the parameter nz
3376: (or the array nnz). By setting these parameters accurately, performance
3377: during matrix assembly can be increased by more than a factor of 50.
3379: Collective on MPI_Comm
3381: Input Parameters:
3382: + comm - MPI communicator, set to PETSC_COMM_SELF
3383: . m - number of rows
3384: . n - number of columns
3385: . nz - number of nonzeros per row (same for all rows)
3386: - nnz - array containing the number of nonzeros in the various rows
3387: (possibly different for each row) or NULL
3389: Output Parameter:
3390: . A - the matrix
3392: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
3393: MatXXXXSetPreallocation() paradgm instead of this routine directly.
3394: [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]
3396: Notes:
3397: If nnz is given then nz is ignored
3399: The AIJ format (also called the Yale sparse matrix format or
3400: compressed row storage), is fully compatible with standard Fortran 77
3401: storage. That is, the stored row and column indices can begin at
3402: either one (as in Fortran) or zero. See the users' manual for details.
3404: Specify the preallocated storage with either nz or nnz (not both).
3405: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
3406: allocation. For large problems you MUST preallocate memory or you
3407: will get TERRIBLE performance, see the users' manual chapter on matrices.
3409: By default, this format uses inodes (identical nodes) when possible, to
3410: improve numerical efficiency of matrix-vector products and solves. We
3411: search for consecutive rows with the same nonzero structure, thereby
3412: reusing matrix information to achieve increased efficiency.
3414: Options Database Keys:
3415: + -mat_no_inode - Do not use inodes
3416: - -mat_inode_limit <limit> - Sets inode limit (max limit=5)
3418: Level: intermediate
3420: .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays()
3422: @*/
3423: PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
3424: {
3428: MatCreate(comm,A);
3429: MatSetSizes(*A,m,n,m,n);
3430: MatSetType(*A,MATSEQAIJ);
3431: MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);
3432: return(0);
3433: }
3437: /*@C
3438: MatSeqAIJSetPreallocation - For good matrix assembly performance
3439: the user should preallocate the matrix storage by setting the parameter nz
3440: (or the array nnz). By setting these parameters accurately, performance
3441: during matrix assembly can be increased by more than a factor of 50.
3443: Collective on MPI_Comm
3445: Input Parameters:
3446: + B - The matrix-free
3447: . nz - number of nonzeros per row (same for all rows)
3448: - nnz - array containing the number of nonzeros in the various rows
3449: (possibly different for each row) or NULL
3451: Notes:
3452: If nnz is given then nz is ignored
3454: The AIJ format (also called the Yale sparse matrix format or
3455: compressed row storage), is fully compatible with standard Fortran 77
3456: storage. That is, the stored row and column indices can begin at
3457: either one (as in Fortran) or zero. See the users' manual for details.
3459: Specify the preallocated storage with either nz or nnz (not both).
3460: Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
3461: allocation. For large problems you MUST preallocate memory or you
3462: will get TERRIBLE performance, see the users' manual chapter on matrices.
3464: You can call MatGetInfo() to get information on how effective the preallocation was;
3465: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
3466: You can also run with the option -info and look for messages with the string
3467: malloc in them to see if additional memory allocation was needed.
3469: Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix
3470: entries or columns indices
3472: By default, this format uses inodes (identical nodes) when possible, to
3473: improve numerical efficiency of matrix-vector products and solves. We
3474: search for consecutive rows with the same nonzero structure, thereby
3475: reusing matrix information to achieve increased efficiency.
3477: Options Database Keys:
3478: + -mat_no_inode - Do not use inodes
3479: . -mat_inode_limit <limit> - Sets inode limit (max limit=5)
3480: - -mat_aij_oneindex - Internally use indexing starting at 1
3481: rather than 0. Note that when calling MatSetValues(),
3482: the user still MUST index entries starting at 0!
3484: Level: intermediate
3486: .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo()
3488: @*/
3489: PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[])
3490: {
3496: PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));
3497: return(0);
3498: }
3502: PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz)
3503: {
3504: Mat_SeqAIJ *b;
3505: PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE;
3507: PetscInt i;
3510: if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
3511: if (nz == MAT_SKIP_ALLOCATION) {
3512: skipallocation = PETSC_TRUE;
3513: nz = 0;
3514: }
3516: PetscLayoutSetUp(B->rmap);
3517: PetscLayoutSetUp(B->cmap);
3519: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
3520: if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz);
3521: if (nnz) {
3522: for (i=0; i<B->rmap->n; i++) {
3523: if (nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]);
3524: if (nnz[i] > B->cmap->n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->cmap->n);
3525: }
3526: }
3528: B->preallocated = PETSC_TRUE;
3530: b = (Mat_SeqAIJ*)B->data;
3532: if (!skipallocation) {
3533: if (!b->imax) {
3534: PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);
3535: PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));
3536: }
3537: if (!nnz) {
3538: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10;
3539: else if (nz < 0) nz = 1;
3540: for (i=0; i<B->rmap->n; i++) b->imax[i] = nz;
3541: nz = nz*B->rmap->n;
3542: } else {
3543: nz = 0;
3544: for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];}
3545: }
3546: /* b->ilen will count nonzeros in each row so far. */
3547: for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0;
3549: /* allocate the matrix space */
3550: MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);
3551: PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);
3552: PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));
3553: b->i[0] = 0;
3554: for (i=1; i<B->rmap->n+1; i++) {
3555: b->i[i] = b->i[i-1] + b->imax[i-1];
3556: }
3557: b->singlemalloc = PETSC_TRUE;
3558: b->free_a = PETSC_TRUE;
3559: b->free_ij = PETSC_TRUE;
3560: #if defined(PETSC_THREADCOMM_ACTIVE)
3561: MatZeroEntries_SeqAIJ(B);
3562: #endif
3563: } else {
3564: b->free_a = PETSC_FALSE;
3565: b->free_ij = PETSC_FALSE;
3566: }
3568: b->nz = 0;
3569: b->maxnz = nz;
3570: B->info.nz_unneeded = (double)b->maxnz;
3571: if (realalloc) {
3572: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
3573: }
3574: return(0);
3575: }
3577: #undef __FUNCT__
3579: /*@
3580: MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format.
3582: Input Parameters:
3583: + B - the matrix
3584: . i - the indices into j for the start of each row (starts with zero)
3585: . j - the column indices for each row (starts with zero) these must be sorted for each row
3586: - v - optional values in the matrix
3588: Level: developer
3590: The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays()
3592: .keywords: matrix, aij, compressed row, sparse, sequential
3594: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ
3595: @*/
3596: PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[])
3597: {
3603: PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));
3604: return(0);
3605: }
3607: #undef __FUNCT__
3609: PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[])
3610: {
3611: PetscInt i;
3612: PetscInt m,n;
3613: PetscInt nz;
3614: PetscInt *nnz, nz_max = 0;
3615: PetscScalar *values;
3619: if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]);
3621: PetscLayoutSetUp(B->rmap);
3622: PetscLayoutSetUp(B->cmap);
3624: MatGetSize(B, &m, &n);
3625: PetscMalloc((m+1) * sizeof(PetscInt), &nnz);
3626: for (i = 0; i < m; i++) {
3627: nz = Ii[i+1]- Ii[i];
3628: nz_max = PetscMax(nz_max, nz);
3629: if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz);
3630: nnz[i] = nz;
3631: }
3632: MatSeqAIJSetPreallocation(B, 0, nnz);
3633: PetscFree(nnz);
3635: if (v) {
3636: values = (PetscScalar*) v;
3637: } else {
3638: PetscMalloc(nz_max*sizeof(PetscScalar), &values);
3639: PetscMemzero(values, nz_max*sizeof(PetscScalar));
3640: }
3642: for (i = 0; i < m; i++) {
3643: nz = Ii[i+1] - Ii[i];
3644: MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);
3645: }
3647: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
3648: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
3650: if (!v) {
3651: PetscFree(values);
3652: }
3653: MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
3654: return(0);
3655: }
3657: #include <../src/mat/impls/dense/seq/dense.h>
3658: #include <petsc-private/kernels/petscaxpy.h>
3662: /*
3663: Computes (B'*A')' since computing B*A directly is untenable
3665: n p p
3666: ( ) ( ) ( )
3667: m ( A ) * n ( B ) = m ( C )
3668: ( ) ( ) ( )
3670: */
3671: PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C)
3672: {
3673: PetscErrorCode ierr;
3674: Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data;
3675: Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data;
3676: Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data;
3677: PetscInt i,n,m,q,p;
3678: const PetscInt *ii,*idx;
3679: const PetscScalar *b,*a,*a_q;
3680: PetscScalar *c,*c_q;
3683: m = A->rmap->n;
3684: n = A->cmap->n;
3685: p = B->cmap->n;
3686: a = sub_a->v;
3687: b = sub_b->a;
3688: c = sub_c->v;
3689: PetscMemzero(c,m*p*sizeof(PetscScalar));
3691: ii = sub_b->i;
3692: idx = sub_b->j;
3693: for (i=0; i<n; i++) {
3694: q = ii[i+1] - ii[i];
3695: while (q-->0) {
3696: c_q = c + m*(*idx);
3697: a_q = a + m*i;
3698: PetscKernelAXPY(c_q,*b,a_q,m);
3699: idx++;
3700: b++;
3701: }
3702: }
3703: return(0);
3704: }
3708: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
3709: {
3711: PetscInt m=A->rmap->n,n=B->cmap->n;
3712: Mat Cmat;
3715: if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n);
3716: MatCreate(PetscObjectComm((PetscObject)A),&Cmat);
3717: MatSetSizes(Cmat,m,n,m,n);
3718: MatSetBlockSizes(Cmat,A->rmap->bs,B->cmap->bs);
3719: MatSetType(Cmat,MATSEQDENSE);
3720: MatSeqDenseSetPreallocation(Cmat,NULL);
3722: Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ;
3724: *C = Cmat;
3725: return(0);
3726: }
3728: /* ----------------------------------------------------------------*/
3731: PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
3732: {
3736: if (scall == MAT_INITIAL_MATRIX) {
3737: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
3738: MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);
3739: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
3740: }
3741: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
3742: MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);
3743: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
3744: return(0);
3745: }
3748: /*MC
3749: MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices,
3750: based on compressed sparse row format.
3752: Options Database Keys:
3753: . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions()
3755: Level: beginner
3757: .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType
3758: M*/
3760: /*MC
3761: MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices.
3763: This matrix type is identical to MATSEQAIJ when constructed with a single process communicator,
3764: and MATMPIAIJ otherwise. As a result, for single process communicators,
3765: MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported
3766: for communicators controlling multiple processes. It is recommended that you call both of
3767: the above preallocation routines for simplicity.
3769: Options Database Keys:
3770: . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions()
3772: Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when
3773: enough exist.
3775: Level: beginner
3777: .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ
3778: M*/
3780: /*MC
3781: MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices.
3783: This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator,
3784: and MATMPIAIJCRL otherwise. As a result, for single process communicators,
3785: MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported
3786: for communicators controlling multiple processes. It is recommended that you call both of
3787: the above preallocation routines for simplicity.
3789: Options Database Keys:
3790: . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions()
3792: Level: beginner
3794: .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL
3795: M*/
3797: #if defined(PETSC_HAVE_PASTIX)
3798: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*);
3799: #endif
3800: #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128)
3801: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat*);
3802: #endif
3803: PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*);
3804: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*);
3805: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*);
3806: extern PetscErrorCode MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscBool*);
3807: #if defined(PETSC_HAVE_MUMPS)
3808: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*);
3809: #endif
3810: #if defined(PETSC_HAVE_SUPERLU)
3811: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*);
3812: #endif
3813: #if defined(PETSC_HAVE_SUPERLU_DIST)
3814: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*);
3815: #endif
3816: #if defined(PETSC_HAVE_UMFPACK)
3817: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*);
3818: #endif
3819: #if defined(PETSC_HAVE_CHOLMOD)
3820: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*);
3821: #endif
3822: #if defined(PETSC_HAVE_LUSOL)
3823: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*);
3824: #endif
3825: #if defined(PETSC_HAVE_MATLAB_ENGINE)
3826: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*);
3827: extern PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*);
3828: extern PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*);
3829: #endif
3830: #if defined(PETSC_HAVE_CLIQUE)
3831: PETSC_EXTERN PetscErrorCode MatGetFactor_aij_clique(Mat,MatFactorType,Mat*);
3832: #endif
3837: /*@C
3838: MatSeqAIJGetArray - gives access to the array where the data for a SeqSeqAIJ matrix is stored
3840: Not Collective
3842: Input Parameter:
3843: . mat - a MATSEQDENSE matrix
3845: Output Parameter:
3846: . array - pointer to the data
3848: Level: intermediate
3850: .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90()
3851: @*/
3852: PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array)
3853: {
3857: PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));
3858: return(0);
3859: }
3863: /*@C
3864: MatSeqAIJRestoreArray - returns access to the array where the data for a SeqSeqAIJ matrix is stored obtained by MatSeqAIJGetArray()
3866: Not Collective
3868: Input Parameters:
3869: . mat - a MATSEQDENSE matrix
3870: . array - pointer to the data
3872: Level: intermediate
3874: .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90()
3875: @*/
3876: PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array)
3877: {
3881: PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));
3882: return(0);
3883: }
3887: PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B)
3888: {
3889: Mat_SeqAIJ *b;
3891: PetscMPIInt size;
3894: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
3895: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1");
3897: PetscNewLog(B,Mat_SeqAIJ,&b);
3899: B->data = (void*)b;
3901: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
3903: b->row = 0;
3904: b->col = 0;
3905: b->icol = 0;
3906: b->reallocs = 0;
3907: b->ignorezeroentries = PETSC_FALSE;
3908: b->roworiented = PETSC_TRUE;
3909: b->nonew = 0;
3910: b->diag = 0;
3911: b->solve_work = 0;
3912: B->spptr = 0;
3913: b->saved_values = 0;
3914: b->idiag = 0;
3915: b->mdiag = 0;
3916: b->ssor_work = 0;
3917: b->omega = 1.0;
3918: b->fshift = 0.0;
3919: b->idiagvalid = PETSC_FALSE;
3920: b->ibdiagvalid = PETSC_FALSE;
3921: b->keepnonzeropattern = PETSC_FALSE;
3922: b->xtoy = 0;
3923: b->XtoY = 0;
3924: B->same_nonzero = PETSC_FALSE;
3926: PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
3927: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);
3928: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);
3930: #if defined(PETSC_HAVE_MATLAB_ENGINE)
3931: PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_matlab_C",MatGetFactor_seqaij_matlab);
3932: PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);
3933: PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);
3934: #endif
3935: #if defined(PETSC_HAVE_PASTIX)
3936: PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_pastix_C",MatGetFactor_seqaij_pastix);
3937: #endif
3938: #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128)
3939: PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_essl_C",MatGetFactor_seqaij_essl);
3940: #endif
3941: #if defined(PETSC_HAVE_SUPERLU)
3942: PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_superlu_C",MatGetFactor_seqaij_superlu);
3943: #endif
3944: #if defined(PETSC_HAVE_SUPERLU_DIST)
3945: PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_superlu_dist_C",MatGetFactor_seqaij_superlu_dist);
3946: #endif
3947: #if defined(PETSC_HAVE_MUMPS)
3948: PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_mumps_C",MatGetFactor_aij_mumps);
3949: #endif
3950: #if defined(PETSC_HAVE_UMFPACK)
3951: PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_umfpack_C",MatGetFactor_seqaij_umfpack);
3952: #endif
3953: #if defined(PETSC_HAVE_CHOLMOD)
3954: PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_cholmod_C",MatGetFactor_seqaij_cholmod);
3955: #endif
3956: #if defined(PETSC_HAVE_LUSOL)
3957: PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_lusol_C",MatGetFactor_seqaij_lusol);
3958: #endif
3959: #if defined(PETSC_HAVE_CLIQUE)
3960: PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_clique_C",MatGetFactor_aij_clique);
3961: #endif
3963: PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_petsc_C",MatGetFactor_seqaij_petsc);
3964: PetscObjectComposeFunction((PetscObject)B,"MatGetFactorAvailable_petsc_C",MatGetFactorAvailable_seqaij_petsc);
3965: PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_bas_C",MatGetFactor_seqaij_bas);
3966: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);
3967: PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);
3968: PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);
3969: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);
3970: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);
3971: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);
3972: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);
3973: PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);
3974: PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);
3975: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);
3976: PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);
3977: PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);
3978: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);
3979: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);
3980: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);
3981: MatCreate_SeqAIJ_Inode(B);
3982: PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
3983: return(0);
3984: }
3988: /*
3989: Given a matrix generated with MatGetFactor() duplicates all the information in A into B
3990: */
3991: PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace)
3992: {
3993: Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data;
3995: PetscInt i,m = A->rmap->n;
3998: c = (Mat_SeqAIJ*)C->data;
4000: C->factortype = A->factortype;
4001: c->row = 0;
4002: c->col = 0;
4003: c->icol = 0;
4004: c->reallocs = 0;
4006: C->assembled = PETSC_TRUE;
4008: PetscLayoutReference(A->rmap,&C->rmap);
4009: PetscLayoutReference(A->cmap,&C->cmap);
4011: PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);
4012: PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));
4013: for (i=0; i<m; i++) {
4014: c->imax[i] = a->imax[i];
4015: c->ilen[i] = a->ilen[i];
4016: }
4018: /* allocate the matrix space */
4019: if (mallocmatspace) {
4020: PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);
4021: PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));
4023: c->singlemalloc = PETSC_TRUE;
4025: PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));
4026: if (m > 0) {
4027: PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));
4028: if (cpvalues == MAT_COPY_VALUES) {
4029: PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));
4030: } else {
4031: PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));
4032: }
4033: }
4034: }
4036: c->ignorezeroentries = a->ignorezeroentries;
4037: c->roworiented = a->roworiented;
4038: c->nonew = a->nonew;
4039: if (a->diag) {
4040: PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);
4041: PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));
4042: for (i=0; i<m; i++) {
4043: c->diag[i] = a->diag[i];
4044: }
4045: } else c->diag = 0;
4047: c->solve_work = 0;
4048: c->saved_values = 0;
4049: c->idiag = 0;
4050: c->ssor_work = 0;
4051: c->keepnonzeropattern = a->keepnonzeropattern;
4052: c->free_a = PETSC_TRUE;
4053: c->free_ij = PETSC_TRUE;
4054: c->xtoy = 0;
4055: c->XtoY = 0;
4057: c->rmax = a->rmax;
4058: c->nz = a->nz;
4059: c->maxnz = a->nz; /* Since we allocate exactly the right amount */
4060: C->preallocated = PETSC_TRUE;
4062: c->compressedrow.use = a->compressedrow.use;
4063: c->compressedrow.nrows = a->compressedrow.nrows;
4064: c->compressedrow.check = a->compressedrow.check;
4065: if (a->compressedrow.use) {
4066: i = a->compressedrow.nrows;
4067: PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);
4068: PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));
4069: PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));
4070: } else {
4071: c->compressedrow.use = PETSC_FALSE;
4072: c->compressedrow.i = NULL;
4073: c->compressedrow.rindex = NULL;
4074: }
4075: C->same_nonzero = A->same_nonzero;
4077: MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);
4078: PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);
4079: return(0);
4080: }
4084: PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
4085: {
4089: MatCreate(PetscObjectComm((PetscObject)A),B);
4090: MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
4091: MatSetBlockSizes(*B,A->rmap->bs,A->cmap->bs);
4092: MatSetType(*B,((PetscObject)A)->type_name);
4093: MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);
4094: return(0);
4095: }
4099: PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer)
4100: {
4101: Mat_SeqAIJ *a;
4103: PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols;
4104: int fd;
4105: PetscMPIInt size;
4106: MPI_Comm comm;
4107: PetscInt bs = 1;
4110: PetscObjectGetComm((PetscObject)viewer,&comm);
4111: MPI_Comm_size(comm,&size);
4112: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor");
4114: PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");
4115: PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);
4116: PetscOptionsEnd();
4117: if (bs > 1) {MatSetBlockSize(newMat,bs);}
4119: PetscViewerBinaryGetDescriptor(viewer,&fd);
4120: PetscBinaryRead(fd,header,4,PETSC_INT);
4121: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file");
4122: M = header[1]; N = header[2]; nz = header[3];
4124: if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ");
4126: /* read in row lengths */
4127: PetscMalloc(M*sizeof(PetscInt),&rowlengths);
4128: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
4130: /* check if sum of rowlengths is same as nz */
4131: for (i=0,sum=0; i< M; i++) sum +=rowlengths[i];
4132: if (sum != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum);
4134: /* set global size if not set already*/
4135: if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) {
4136: MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);
4137: } else {
4138: /* if sizes and type are already set, check if the vector global sizes are correct */
4139: MatGetSize(newMat,&rows,&cols);
4140: if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */
4141: MatGetLocalSize(newMat,&rows,&cols);
4142: }
4143: if (M != rows || N != cols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols);
4144: }
4145: MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);
4146: a = (Mat_SeqAIJ*)newMat->data;
4148: PetscBinaryRead(fd,a->j,nz,PETSC_INT);
4150: /* read in nonzero values */
4151: PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);
4153: /* set matrix "i" values */
4154: a->i[0] = 0;
4155: for (i=1; i<= M; i++) {
4156: a->i[i] = a->i[i-1] + rowlengths[i-1];
4157: a->ilen[i-1] = rowlengths[i-1];
4158: }
4159: PetscFree(rowlengths);
4161: MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);
4162: MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);
4163: return(0);
4164: }
4168: PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg)
4169: {
4170: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data;
4172: #if defined(PETSC_USE_COMPLEX)
4173: PetscInt k;
4174: #endif
4177: /* If the matrix dimensions are not equal,or no of nonzeros */
4178: if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) {
4179: *flg = PETSC_FALSE;
4180: return(0);
4181: }
4183: /* if the a->i are the same */
4184: PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);
4185: if (!*flg) return(0);
4187: /* if a->j are the same */
4188: PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);
4189: if (!*flg) return(0);
4191: /* if a->a are the same */
4192: #if defined(PETSC_USE_COMPLEX)
4193: for (k=0; k<a->nz; k++) {
4194: if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) {
4195: *flg = PETSC_FALSE;
4196: return(0);
4197: }
4198: }
4199: #else
4200: PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);
4201: #endif
4202: return(0);
4203: }
4207: /*@
4208: MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format)
4209: provided by the user.
4211: Collective on MPI_Comm
4213: Input Parameters:
4214: + comm - must be an MPI communicator of size 1
4215: . m - number of rows
4216: . n - number of columns
4217: . i - row indices
4218: . j - column indices
4219: - a - matrix values
4221: Output Parameter:
4222: . mat - the matrix
4224: Level: intermediate
4226: Notes:
4227: The i, j, and a arrays are not copied by this routine, the user must free these arrays
4228: once the matrix is destroyed and not before
4230: You cannot set new nonzero locations into this matrix, that will generate an error.
4232: The i and j indices are 0 based
4234: The format which is used for the sparse matrix input, is equivalent to a
4235: row-major ordering.. i.e for the following matrix, the input data expected is
4236: as shown:
4238: 1 0 0
4239: 2 0 3
4240: 4 5 6
4242: i = {0,1,3,6} [size = nrow+1 = 3+1]
4243: j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row
4244: v = {1,2,3,4,5,6} [size = nz = 6]
4247: .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR()
4249: @*/
4250: PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat)
4251: {
4253: PetscInt ii;
4254: Mat_SeqAIJ *aij;
4255: #if defined(PETSC_USE_DEBUG)
4256: PetscInt jj;
4257: #endif
4260: if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0");
4261: MatCreate(comm,mat);
4262: MatSetSizes(*mat,m,n,m,n);
4263: /* MatSetBlockSizes(*mat,,); */
4264: MatSetType(*mat,MATSEQAIJ);
4265: MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);
4266: aij = (Mat_SeqAIJ*)(*mat)->data;
4267: PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);
4269: aij->i = i;
4270: aij->j = j;
4271: aij->a = a;
4272: aij->singlemalloc = PETSC_FALSE;
4273: aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
4274: aij->free_a = PETSC_FALSE;
4275: aij->free_ij = PETSC_FALSE;
4277: for (ii=0; ii<m; ii++) {
4278: aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii];
4279: #if defined(PETSC_USE_DEBUG)
4280: if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]);
4281: for (jj=i[ii]+1; jj<i[ii+1]; jj++) {
4282: if (j[jj] < j[jj-1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is not sorted",jj-i[ii],j[jj],ii);
4283: if (j[jj] == j[jj]-1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is identical to previous entry",jj-i[ii],j[jj],ii);
4284: }
4285: #endif
4286: }
4287: #if defined(PETSC_USE_DEBUG)
4288: for (ii=0; ii<aij->i[m]; ii++) {
4289: if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]);
4290: if (j[ii] > n - 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]);
4291: }
4292: #endif
4294: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
4295: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
4296: return(0);
4297: }
4300: /*@C
4301: MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format)
4302: provided by the user.
4304: Collective on MPI_Comm
4306: Input Parameters:
4307: + comm - must be an MPI communicator of size 1
4308: . m - number of rows
4309: . n - number of columns
4310: . i - row indices
4311: . j - column indices
4312: . a - matrix values
4313: . nz - number of nonzeros
4314: - idx - 0 or 1 based
4316: Output Parameter:
4317: . mat - the matrix
4319: Level: intermediate
4321: Notes:
4322: The i and j indices are 0 based
4324: The format which is used for the sparse matrix input, is equivalent to a
4325: row-major ordering.. i.e for the following matrix, the input data expected is
4326: as shown:
4328: 1 0 0
4329: 2 0 3
4330: 4 5 6
4332: i = {0,1,1,2,2,2}
4333: j = {0,0,2,0,1,2}
4334: v = {1,2,3,4,5,6}
4337: .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR()
4339: @*/
4340: PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat,PetscInt nz,PetscBool idx)
4341: {
4343: PetscInt ii, *nnz, one = 1,row,col;
4347: PetscMalloc(m*sizeof(PetscInt),&nnz);
4348: PetscMemzero(nnz,m*sizeof(PetscInt));
4349: for (ii = 0; ii < nz; ii++) {
4350: nnz[i[ii] - !!idx] += 1;
4351: }
4352: MatCreate(comm,mat);
4353: MatSetSizes(*mat,m,n,m,n);
4354: MatSetType(*mat,MATSEQAIJ);
4355: MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);
4356: for (ii = 0; ii < nz; ii++) {
4357: if (idx) {
4358: row = i[ii] - 1;
4359: col = j[ii] - 1;
4360: } else {
4361: row = i[ii];
4362: col = j[ii];
4363: }
4364: MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);
4365: }
4366: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
4367: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
4368: PetscFree(nnz);
4369: return(0);
4370: }
4374: PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring)
4375: {
4377: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
4380: if (coloring->ctype == IS_COLORING_GLOBAL) {
4381: ISColoringReference(coloring);
4382: a->coloring = coloring;
4383: } else if (coloring->ctype == IS_COLORING_GHOSTED) {
4384: PetscInt i,*larray;
4385: ISColoring ocoloring;
4386: ISColoringValue *colors;
4388: /* set coloring for diagonal portion */
4389: PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);
4390: for (i=0; i<A->cmap->n; i++) larray[i] = i;
4391: ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,NULL,larray);
4392: PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);
4393: for (i=0; i<A->cmap->n; i++) colors[i] = coloring->colors[larray[i]];
4394: PetscFree(larray);
4395: ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);
4396: a->coloring = ocoloring;
4397: }
4398: return(0);
4399: }
4403: PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues)
4404: {
4405: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
4406: PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j;
4407: MatScalar *v = a->a;
4408: PetscScalar *values = (PetscScalar*)advalues;
4409: ISColoringValue *color;
4412: if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix");
4413: color = a->coloring->colors;
4414: /* loop over rows */
4415: for (i=0; i<m; i++) {
4416: nz = ii[i+1] - ii[i];
4417: /* loop over columns putting computed value into matrix */
4418: for (j=0; j<nz; j++) *v++ = values[color[*jj++]];
4419: values += nl; /* jump to next row of derivatives */
4420: }
4421: return(0);
4422: }
4426: PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A)
4427: {
4428: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
4432: a->idiagvalid = PETSC_FALSE;
4433: a->ibdiagvalid = PETSC_FALSE;
4435: MatSeqAIJInvalidateDiagonal_Inode(A);
4436: return(0);
4437: }
4439: /*
4440: Special version for direct calls from Fortran
4441: */
4442: #include <petsc-private/fortranimpl.h>
4443: #if defined(PETSC_HAVE_FORTRAN_CAPS)
4444: #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ
4445: #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
4446: #define matsetvaluesseqaij_ matsetvaluesseqaij
4447: #endif
4449: /* Change these macros so can be used in void function */
4450: #undef CHKERRQ
4451: #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr)
4452: #undef SETERRQ2
4453: #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr)
4454: #undef SETERRQ3
4455: #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr)
4459: PETSC_EXTERN void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr)
4460: {
4461: Mat A = *AA;
4462: PetscInt m = *mm, n = *nn;
4463: InsertMode is = *isis;
4464: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
4465: PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N;
4466: PetscInt *imax,*ai,*ailen;
4468: PetscInt *aj,nonew = a->nonew,lastcol = -1;
4469: MatScalar *ap,value,*aa;
4470: PetscBool ignorezeroentries = a->ignorezeroentries;
4471: PetscBool roworiented = a->roworiented;
4474: MatCheckPreallocated(A,1);
4475: imax = a->imax;
4476: ai = a->i;
4477: ailen = a->ilen;
4478: aj = a->j;
4479: aa = a->a;
4481: for (k=0; k<m; k++) { /* loop over added rows */
4482: row = im[k];
4483: if (row < 0) continue;
4484: #if defined(PETSC_USE_DEBUG)
4485: if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large");
4486: #endif
4487: rp = aj + ai[row]; ap = aa + ai[row];
4488: rmax = imax[row]; nrow = ailen[row];
4489: low = 0;
4490: high = nrow;
4491: for (l=0; l<n; l++) { /* loop over added columns */
4492: if (in[l] < 0) continue;
4493: #if defined(PETSC_USE_DEBUG)
4494: if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large");
4495: #endif
4496: col = in[l];
4497: if (roworiented) value = v[l + k*n];
4498: else value = v[k + l*m];
4500: if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue;
4502: if (col <= lastcol) low = 0;
4503: else high = nrow;
4504: lastcol = col;
4505: while (high-low > 5) {
4506: t = (low+high)/2;
4507: if (rp[t] > col) high = t;
4508: else low = t;
4509: }
4510: for (i=low; i<high; i++) {
4511: if (rp[i] > col) break;
4512: if (rp[i] == col) {
4513: if (is == ADD_VALUES) ap[i] += value;
4514: else ap[i] = value;
4515: goto noinsert;
4516: }
4517: }
4518: if (value == 0.0 && ignorezeroentries) goto noinsert;
4519: if (nonew == 1) goto noinsert;
4520: if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix");
4521: MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
4522: N = nrow++ - 1; a->nz++; high++;
4523: /* shift up all the later entries in this row */
4524: for (ii=N; ii>=i; ii--) {
4525: rp[ii+1] = rp[ii];
4526: ap[ii+1] = ap[ii];
4527: }
4528: rp[i] = col;
4529: ap[i] = value;
4530: noinsert:;
4531: low = i + 1;
4532: }
4533: ailen[row] = nrow;
4534: }
4535: A->same_nonzero = PETSC_FALSE;
4536: PetscFunctionReturnVoid();
4537: }