Actual source code: dense.c
petsc-3.8.4 2018-03-24
2: /*
3: Defines the basic matrix operations for sequential dense.
4: */
6: #include <../src/mat/impls/dense/seq/dense.h>
7: #include <petscblaslapack.h>
9: #include <../src/mat/impls/aij/seq/aij.h>
11: static PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian)
12: {
13: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
14: PetscInt j, k, n = A->rmap->n;
17: if (A->rmap->n != A->cmap->n) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot symmetrize a rectangular matrix");
18: if (!hermitian) {
19: for (k=0;k<n;k++) {
20: for (j=k;j<n;j++) {
21: mat->v[j*mat->lda + k] = mat->v[k*mat->lda + j];
22: }
23: }
24: } else {
25: for (k=0;k<n;k++) {
26: for (j=k;j<n;j++) {
27: mat->v[j*mat->lda + k] = PetscConj(mat->v[k*mat->lda + j]);
28: }
29: }
30: }
31: return(0);
32: }
34: PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A)
35: {
36: #if defined(PETSC_MISSING_LAPACK_POTRF)
38: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
39: #else
40: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
42: PetscBLASInt info,n;
45: if (!A->rmap->n || !A->cmap->n) return(0);
46: PetscBLASIntCast(A->cmap->n,&n);
47: if (A->factortype == MAT_FACTOR_LU) {
48: if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
49: if (!mat->fwork) {
50: mat->lfwork = n;
51: PetscMalloc1(mat->lfwork,&mat->fwork);
52: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
53: }
54: PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
55: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0); /* TODO CHECK FLOPS */
56: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
57: if (A->spd) {
58: PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info));
59: MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
60: #if defined (PETSC_USE_COMPLEX)
61: } else if (A->hermitian) {
62: if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
63: if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present");
64: PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
65: MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
66: #endif
67: } else { /* symmetric case */
68: if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
69: if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present");
70: PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
71: MatSeqDenseSymmetrize_Private(A,PETSC_FALSE);
72: }
73: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad Inversion: zero pivot in row %D",(PetscInt)info-1);
74: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0); /* TODO CHECK FLOPS */
75: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
76: #endif
78: A->ops->solve = NULL;
79: A->ops->matsolve = NULL;
80: A->ops->solvetranspose = NULL;
81: A->ops->matsolvetranspose = NULL;
82: A->ops->solveadd = NULL;
83: A->ops->solvetransposeadd = NULL;
84: A->factortype = MAT_FACTOR_NONE;
85: PetscFree(A->solvertype);
86: return(0);
87: }
89: PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
90: {
91: PetscErrorCode ierr;
92: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
93: PetscInt m = l->lda, n = A->cmap->n,r = A->rmap->n, i,j;
94: PetscScalar *slot,*bb;
95: const PetscScalar *xx;
98: #if defined(PETSC_USE_DEBUG)
99: for (i=0; i<N; i++) {
100: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
101: if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n);
102: if (rows[i] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Col %D requested to be zeroed greater than or equal number of cols %D",rows[i],A->cmap->n);
103: }
104: #endif
106: /* fix right hand side if needed */
107: if (x && b) {
108: VecGetArrayRead(x,&xx);
109: VecGetArray(b,&bb);
110: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
111: VecRestoreArrayRead(x,&xx);
112: VecRestoreArray(b,&bb);
113: }
115: for (i=0; i<N; i++) {
116: slot = l->v + rows[i]*m;
117: PetscMemzero(slot,r*sizeof(PetscScalar));
118: }
119: for (i=0; i<N; i++) {
120: slot = l->v + rows[i];
121: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
122: }
123: if (diag != 0.0) {
124: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
125: for (i=0; i<N; i++) {
126: slot = l->v + (m+1)*rows[i];
127: *slot = diag;
128: }
129: }
130: return(0);
131: }
133: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
134: {
135: Mat_SeqDense *c = (Mat_SeqDense*)(C->data);
139: MatMatMultNumeric_SeqDense_SeqDense(A,P,c->ptapwork);
140: MatTransposeMatMultNumeric_SeqDense_SeqDense(P,c->ptapwork,C);
141: return(0);
142: }
144: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat *C)
145: {
146: Mat_SeqDense *c;
150: MatCreateSeqDense(PetscObjectComm((PetscObject)A),P->cmap->N,P->cmap->N,NULL,C);
151: c = (Mat_SeqDense*)((*C)->data);
152: MatCreateSeqDense(PetscObjectComm((PetscObject)A),A->rmap->N,P->cmap->N,NULL,&c->ptapwork);
153: return(0);
154: }
156: PETSC_INTERN PetscErrorCode MatPtAP_SeqDense_SeqDense(Mat A,Mat P,MatReuse reuse,PetscReal fill,Mat *C)
157: {
161: if (reuse == MAT_INITIAL_MATRIX) {
162: MatPtAPSymbolic_SeqDense_SeqDense(A,P,fill,C);
163: }
164: PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
165: (*(*C)->ops->ptapnumeric)(A,P,*C);
166: PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
167: return(0);
168: }
170: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
171: {
172: Mat B = NULL;
173: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
174: Mat_SeqDense *b;
176: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
177: MatScalar *av=a->a;
178: PetscBool isseqdense;
181: if (reuse == MAT_REUSE_MATRIX) {
182: PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);
183: if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type);
184: }
185: if (reuse != MAT_REUSE_MATRIX) {
186: MatCreate(PetscObjectComm((PetscObject)A),&B);
187: MatSetSizes(B,m,n,m,n);
188: MatSetType(B,MATSEQDENSE);
189: MatSeqDenseSetPreallocation(B,NULL);
190: b = (Mat_SeqDense*)(B->data);
191: } else {
192: b = (Mat_SeqDense*)((*newmat)->data);
193: PetscMemzero(b->v,m*n*sizeof(PetscScalar));
194: }
195: for (i=0; i<m; i++) {
196: PetscInt j;
197: for (j=0;j<ai[1]-ai[0];j++) {
198: b->v[*aj*m+i] = *av;
199: aj++;
200: av++;
201: }
202: ai++;
203: }
205: if (reuse == MAT_INPLACE_MATRIX) {
206: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
207: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
208: MatHeaderReplace(A,&B);
209: } else {
210: if (B) *newmat = B;
211: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
212: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
213: }
214: return(0);
215: }
217: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
218: {
219: Mat B;
220: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
222: PetscInt i, j;
223: PetscInt *rows, *nnz;
224: MatScalar *aa = a->v, *vals;
227: MatCreate(PetscObjectComm((PetscObject)A),&B);
228: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
229: MatSetType(B,MATSEQAIJ);
230: PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
231: for (j=0; j<A->cmap->n; j++) {
232: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
233: aa += a->lda;
234: }
235: MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
236: aa = a->v;
237: for (j=0; j<A->cmap->n; j++) {
238: PetscInt numRows = 0;
239: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
240: MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
241: aa += a->lda;
242: }
243: PetscFree3(rows,nnz,vals);
244: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
245: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
247: if (reuse == MAT_INPLACE_MATRIX) {
248: MatHeaderReplace(A,&B);
249: } else {
250: *newmat = B;
251: }
252: return(0);
253: }
255: static PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
256: {
257: Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
258: PetscScalar oalpha = alpha;
259: PetscInt j;
260: PetscBLASInt N,m,ldax,lday,one = 1;
264: PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
265: PetscBLASIntCast(X->rmap->n,&m);
266: PetscBLASIntCast(x->lda,&ldax);
267: PetscBLASIntCast(y->lda,&lday);
268: if (ldax>m || lday>m) {
269: for (j=0; j<X->cmap->n; j++) {
270: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one));
271: }
272: } else {
273: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one));
274: }
275: PetscObjectStateIncrease((PetscObject)Y);
276: PetscLogFlops(PetscMax(2*N-1,0));
277: return(0);
278: }
280: static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
281: {
282: PetscInt N = A->rmap->n*A->cmap->n;
285: info->block_size = 1.0;
286: info->nz_allocated = (double)N;
287: info->nz_used = (double)N;
288: info->nz_unneeded = (double)0;
289: info->assemblies = (double)A->num_ass;
290: info->mallocs = 0;
291: info->memory = ((PetscObject)A)->mem;
292: info->fill_ratio_given = 0;
293: info->fill_ratio_needed = 0;
294: info->factor_mallocs = 0;
295: return(0);
296: }
298: static PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
299: {
300: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
301: PetscScalar oalpha = alpha;
303: PetscBLASInt one = 1,j,nz,lda;
306: PetscBLASIntCast(a->lda,&lda);
307: if (lda>A->rmap->n) {
308: PetscBLASIntCast(A->rmap->n,&nz);
309: for (j=0; j<A->cmap->n; j++) {
310: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
311: }
312: } else {
313: PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
314: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
315: }
316: PetscLogFlops(nz);
317: return(0);
318: }
320: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl)
321: {
322: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
323: PetscInt i,j,m = A->rmap->n,N;
324: PetscScalar *v = a->v;
327: *fl = PETSC_FALSE;
328: if (A->rmap->n != A->cmap->n) return(0);
329: N = a->lda;
331: for (i=0; i<m; i++) {
332: for (j=i+1; j<m; j++) {
333: if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
334: }
335: }
336: *fl = PETSC_TRUE;
337: return(0);
338: }
340: static PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
341: {
342: Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l;
344: PetscInt lda = (PetscInt)mat->lda,j,m;
347: PetscLayoutReference(A->rmap,&newi->rmap);
348: PetscLayoutReference(A->cmap,&newi->cmap);
349: MatSeqDenseSetPreallocation(newi,NULL);
350: if (cpvalues == MAT_COPY_VALUES) {
351: l = (Mat_SeqDense*)newi->data;
352: if (lda>A->rmap->n) {
353: m = A->rmap->n;
354: for (j=0; j<A->cmap->n; j++) {
355: PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
356: }
357: } else {
358: PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
359: }
360: }
361: newi->assembled = PETSC_TRUE;
362: return(0);
363: }
365: static PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
366: {
370: MatCreate(PetscObjectComm((PetscObject)A),newmat);
371: MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
372: MatSetType(*newmat,((PetscObject)A)->type_name);
373: MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
374: return(0);
375: }
378: static PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*);
380: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
381: {
382: MatFactorInfo info;
386: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
387: MatLUFactor_SeqDense(fact,0,0,&info);
388: return(0);
389: }
391: static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
392: {
393: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
394: PetscErrorCode ierr;
395: const PetscScalar *x;
396: PetscScalar *y;
397: PetscBLASInt one = 1,info,m;
400: PetscBLASIntCast(A->rmap->n,&m);
401: VecGetArrayRead(xx,&x);
402: VecGetArray(yy,&y);
403: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
404: if (A->factortype == MAT_FACTOR_LU) {
405: #if defined(PETSC_MISSING_LAPACK_GETRS)
406: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
407: #else
408: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
409: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
410: #endif
411: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
412: #if defined(PETSC_MISSING_LAPACK_POTRS)
413: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
414: #else
415: if (A->spd) {
416: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
417: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
418: #if defined (PETSC_USE_COMPLEX)
419: } else if (A->hermitian) {
420: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
421: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
422: #endif
423: } else { /* symmetric case */
424: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
425: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
426: }
427: #endif
428: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
429: VecRestoreArrayRead(xx,&x);
430: VecRestoreArray(yy,&y);
431: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
432: return(0);
433: }
435: static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
436: {
437: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
439: PetscScalar *b,*x;
440: PetscInt n;
441: PetscBLASInt nrhs,info,m;
442: PetscBool flg;
445: PetscBLASIntCast(A->rmap->n,&m);
446: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
447: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
448: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
449: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
451: MatGetSize(B,NULL,&n);
452: PetscBLASIntCast(n,&nrhs);
453: MatDenseGetArray(B,&b);
454: MatDenseGetArray(X,&x);
456: PetscMemcpy(x,b,m*nrhs*sizeof(PetscScalar));
458: if (A->factortype == MAT_FACTOR_LU) {
459: #if defined(PETSC_MISSING_LAPACK_GETRS)
460: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
461: #else
462: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
463: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
464: #endif
465: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
466: #if defined(PETSC_MISSING_LAPACK_POTRS)
467: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
468: #else
469: if (A->spd) {
470: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
471: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
472: #if defined (PETSC_USE_COMPLEX)
473: } else if (A->hermitian) {
474: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
475: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
476: #endif
477: } else { /* symmetric case */
478: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
479: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
480: }
481: #endif
482: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
484: MatDenseRestoreArray(B,&b);
485: MatDenseRestoreArray(X,&x);
486: PetscLogFlops(nrhs*(2.0*m*m - m));
487: return(0);
488: }
490: static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
491: {
492: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
493: PetscErrorCode ierr;
494: const PetscScalar *x;
495: PetscScalar *y;
496: PetscBLASInt one = 1,info,m;
499: PetscBLASIntCast(A->rmap->n,&m);
500: VecGetArrayRead(xx,&x);
501: VecGetArray(yy,&y);
502: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
503: if (A->factortype == MAT_FACTOR_LU) {
504: #if defined(PETSC_MISSING_LAPACK_GETRS)
505: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
506: #else
507: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
508: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
509: #endif
510: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
511: #if defined(PETSC_MISSING_LAPACK_POTRS)
512: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
513: #else
514: if (A->spd) {
515: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
516: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
517: #if defined (PETSC_USE_COMPLEX)
518: } else if (A->hermitian) {
519: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatSolveTranspose with Cholesky/Hermitian not available");
520: #endif
521: } else { /* symmetric case */
522: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
523: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
524: }
525: #endif
526: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
527: VecRestoreArrayRead(xx,&x);
528: VecRestoreArray(yy,&y);
529: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
530: return(0);
531: }
533: static PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
534: {
535: PetscErrorCode ierr;
536: const PetscScalar *x;
537: PetscScalar *y,sone = 1.0;
538: Vec tmp = 0;
541: VecGetArrayRead(xx,&x);
542: VecGetArray(yy,&y);
543: if (!A->rmap->n || !A->cmap->n) return(0);
544: if (yy == zz) {
545: VecDuplicate(yy,&tmp);
546: PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
547: VecCopy(yy,tmp);
548: }
549: MatSolve_SeqDense(A,xx,yy);
550: if (tmp) {
551: VecAXPY(yy,sone,tmp);
552: VecDestroy(&tmp);
553: } else {
554: VecAXPY(yy,sone,zz);
555: }
556: VecRestoreArrayRead(xx,&x);
557: VecRestoreArray(yy,&y);
558: PetscLogFlops(A->cmap->n);
559: return(0);
560: }
562: static PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
563: {
564: PetscErrorCode ierr;
565: const PetscScalar *x;
566: PetscScalar *y,sone = 1.0;
567: Vec tmp = NULL;
570: if (!A->rmap->n || !A->cmap->n) return(0);
571: VecGetArrayRead(xx,&x);
572: VecGetArray(yy,&y);
573: if (yy == zz) {
574: VecDuplicate(yy,&tmp);
575: PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
576: VecCopy(yy,tmp);
577: }
578: MatSolveTranspose_SeqDense(A,xx,yy);
579: if (tmp) {
580: VecAXPY(yy,sone,tmp);
581: VecDestroy(&tmp);
582: } else {
583: VecAXPY(yy,sone,zz);
584: }
585: VecRestoreArrayRead(xx,&x);
586: VecRestoreArray(yy,&y);
587: PetscLogFlops(A->cmap->n);
588: return(0);
589: }
591: /* ---------------------------------------------------------------*/
592: /* COMMENT: I have chosen to hide row permutation in the pivots,
593: rather than put it in the Mat->row slot.*/
594: static PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
595: {
596: #if defined(PETSC_MISSING_LAPACK_GETRF)
598: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
599: #else
600: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
602: PetscBLASInt n,m,info;
605: PetscBLASIntCast(A->cmap->n,&n);
606: PetscBLASIntCast(A->rmap->n,&m);
607: if (!mat->pivots) {
608: PetscMalloc1(A->rmap->n,&mat->pivots);
609: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
610: }
611: if (!A->rmap->n || !A->cmap->n) return(0);
612: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
613: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
614: PetscFPTrapPop();
616: if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
617: if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
619: A->ops->solve = MatSolve_SeqDense;
620: A->ops->matsolve = MatMatSolve_SeqDense;
621: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
622: A->ops->solveadd = MatSolveAdd_SeqDense;
623: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
624: A->factortype = MAT_FACTOR_LU;
626: PetscFree(A->solvertype);
627: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
629: PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
630: #endif
631: return(0);
632: }
634: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
635: static PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
636: {
637: #if defined(PETSC_MISSING_LAPACK_POTRF)
639: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
640: #else
641: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
643: PetscBLASInt info,n;
646: PetscBLASIntCast(A->cmap->n,&n);
647: if (!A->rmap->n || !A->cmap->n) return(0);
648: if (A->spd) {
649: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
650: #if defined (PETSC_USE_COMPLEX)
651: } else if (A->hermitian) {
652: if (!mat->pivots) {
653: PetscMalloc1(A->rmap->n,&mat->pivots);
654: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
655: }
656: if (!mat->fwork) {
657: PetscScalar dummy;
659: mat->lfwork = -1;
660: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
661: mat->lfwork = (PetscInt)PetscRealPart(dummy);
662: PetscMalloc1(mat->lfwork,&mat->fwork);
663: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
664: }
665: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
666: #endif
667: } else { /* symmetric case */
668: if (!mat->pivots) {
669: PetscMalloc1(A->rmap->n,&mat->pivots);
670: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
671: }
672: if (!mat->fwork) {
673: PetscScalar dummy;
675: mat->lfwork = -1;
676: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
677: mat->lfwork = (PetscInt)PetscRealPart(dummy);
678: PetscMalloc1(mat->lfwork,&mat->fwork);
679: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
680: }
681: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
682: }
683: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
685: A->ops->solve = MatSolve_SeqDense;
686: A->ops->matsolve = MatMatSolve_SeqDense;
687: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
688: A->ops->solveadd = MatSolveAdd_SeqDense;
689: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
690: A->factortype = MAT_FACTOR_CHOLESKY;
692: PetscFree(A->solvertype);
693: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
695: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
696: #endif
697: return(0);
698: }
701: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
702: {
704: MatFactorInfo info;
707: info.fill = 1.0;
709: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
710: MatCholeskyFactor_SeqDense(fact,0,&info);
711: return(0);
712: }
714: static PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
715: {
717: fact->assembled = PETSC_TRUE;
718: fact->preallocated = PETSC_TRUE;
719: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
720: fact->ops->solve = MatSolve_SeqDense;
721: fact->ops->matsolve = MatMatSolve_SeqDense;
722: fact->ops->solvetranspose = MatSolveTranspose_SeqDense;
723: fact->ops->solveadd = MatSolveAdd_SeqDense;
724: fact->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
725: return(0);
726: }
728: static PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
729: {
731: fact->preallocated = PETSC_TRUE;
732: fact->assembled = PETSC_TRUE;
733: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
734: fact->ops->solve = MatSolve_SeqDense;
735: fact->ops->matsolve = MatMatSolve_SeqDense;
736: fact->ops->solvetranspose = MatSolveTranspose_SeqDense;
737: fact->ops->solveadd = MatSolveAdd_SeqDense;
738: fact->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
739: return(0);
740: }
742: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
743: {
747: MatCreate(PetscObjectComm((PetscObject)A),fact);
748: MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
749: MatSetType(*fact,((PetscObject)A)->type_name);
750: if (ftype == MAT_FACTOR_LU) {
751: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
752: } else {
753: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
754: }
755: (*fact)->factortype = ftype;
757: PetscFree((*fact)->solvertype);
758: PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
759: return(0);
760: }
762: /* ------------------------------------------------------------------*/
763: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
764: {
765: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
766: PetscScalar *x,*v = mat->v,zero = 0.0,xt;
767: const PetscScalar *b;
768: PetscErrorCode ierr;
769: PetscInt m = A->rmap->n,i;
770: PetscBLASInt o = 1,bm;
773: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
774: PetscBLASIntCast(m,&bm);
775: if (flag & SOR_ZERO_INITIAL_GUESS) {
776: /* this is a hack fix, should have another version without the second BLASdotu */
777: VecSet(xx,zero);
778: }
779: VecGetArray(xx,&x);
780: VecGetArrayRead(bb,&b);
781: its = its*lits;
782: if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
783: while (its--) {
784: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
785: for (i=0; i<m; i++) {
786: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
787: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
788: }
789: }
790: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
791: for (i=m-1; i>=0; i--) {
792: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
793: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
794: }
795: }
796: }
797: VecRestoreArrayRead(bb,&b);
798: VecRestoreArray(xx,&x);
799: return(0);
800: }
802: /* -----------------------------------------------------------------*/
803: PETSC_INTERN PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
804: {
805: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
806: const PetscScalar *v = mat->v,*x;
807: PetscScalar *y;
808: PetscErrorCode ierr;
809: PetscBLASInt m, n,_One=1;
810: PetscScalar _DOne=1.0,_DZero=0.0;
813: PetscBLASIntCast(A->rmap->n,&m);
814: PetscBLASIntCast(A->cmap->n,&n);
815: VecGetArrayRead(xx,&x);
816: VecGetArray(yy,&y);
817: if (!A->rmap->n || !A->cmap->n) {
818: PetscBLASInt i;
819: for (i=0; i<n; i++) y[i] = 0.0;
820: } else {
821: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
822: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
823: }
824: VecRestoreArrayRead(xx,&x);
825: VecRestoreArray(yy,&y);
826: return(0);
827: }
829: PETSC_INTERN PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
830: {
831: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
832: PetscScalar *y,_DOne=1.0,_DZero=0.0;
833: PetscErrorCode ierr;
834: PetscBLASInt m, n, _One=1;
835: const PetscScalar *v = mat->v,*x;
838: PetscBLASIntCast(A->rmap->n,&m);
839: PetscBLASIntCast(A->cmap->n,&n);
840: VecGetArrayRead(xx,&x);
841: VecGetArray(yy,&y);
842: if (!A->rmap->n || !A->cmap->n) {
843: PetscBLASInt i;
844: for (i=0; i<m; i++) y[i] = 0.0;
845: } else {
846: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
847: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
848: }
849: VecRestoreArrayRead(xx,&x);
850: VecRestoreArray(yy,&y);
851: return(0);
852: }
854: PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
855: {
856: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
857: const PetscScalar *v = mat->v,*x;
858: PetscScalar *y,_DOne=1.0;
859: PetscErrorCode ierr;
860: PetscBLASInt m, n, _One=1;
863: PetscBLASIntCast(A->rmap->n,&m);
864: PetscBLASIntCast(A->cmap->n,&n);
865: if (!A->rmap->n || !A->cmap->n) return(0);
866: if (zz != yy) {VecCopy(zz,yy);}
867: VecGetArrayRead(xx,&x);
868: VecGetArray(yy,&y);
869: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
870: VecRestoreArrayRead(xx,&x);
871: VecRestoreArray(yy,&y);
872: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
873: return(0);
874: }
876: static PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
877: {
878: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
879: const PetscScalar *v = mat->v,*x;
880: PetscScalar *y;
881: PetscErrorCode ierr;
882: PetscBLASInt m, n, _One=1;
883: PetscScalar _DOne=1.0;
886: PetscBLASIntCast(A->rmap->n,&m);
887: PetscBLASIntCast(A->cmap->n,&n);
888: if (!A->rmap->n || !A->cmap->n) return(0);
889: if (zz != yy) {VecCopy(zz,yy);}
890: VecGetArrayRead(xx,&x);
891: VecGetArray(yy,&y);
892: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
893: VecRestoreArrayRead(xx,&x);
894: VecRestoreArray(yy,&y);
895: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
896: return(0);
897: }
899: /* -----------------------------------------------------------------*/
900: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
901: {
902: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
903: PetscScalar *v;
905: PetscInt i;
908: *ncols = A->cmap->n;
909: if (cols) {
910: PetscMalloc1(A->cmap->n+1,cols);
911: for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
912: }
913: if (vals) {
914: PetscMalloc1(A->cmap->n+1,vals);
915: v = mat->v + row;
916: for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
917: }
918: return(0);
919: }
921: static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
922: {
926: if (cols) {PetscFree(*cols);}
927: if (vals) {PetscFree(*vals); }
928: return(0);
929: }
930: /* ----------------------------------------------------------------*/
931: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
932: {
933: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
934: PetscInt i,j,idx=0;
937: if (!mat->roworiented) {
938: if (addv == INSERT_VALUES) {
939: for (j=0; j<n; j++) {
940: if (indexn[j] < 0) {idx += m; continue;}
941: #if defined(PETSC_USE_DEBUG)
942: if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
943: #endif
944: for (i=0; i<m; i++) {
945: if (indexm[i] < 0) {idx++; continue;}
946: #if defined(PETSC_USE_DEBUG)
947: if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
948: #endif
949: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
950: }
951: }
952: } else {
953: for (j=0; j<n; j++) {
954: if (indexn[j] < 0) {idx += m; continue;}
955: #if defined(PETSC_USE_DEBUG)
956: if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
957: #endif
958: for (i=0; i<m; i++) {
959: if (indexm[i] < 0) {idx++; continue;}
960: #if defined(PETSC_USE_DEBUG)
961: if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
962: #endif
963: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
964: }
965: }
966: }
967: } else {
968: if (addv == INSERT_VALUES) {
969: for (i=0; i<m; i++) {
970: if (indexm[i] < 0) { idx += n; continue;}
971: #if defined(PETSC_USE_DEBUG)
972: if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
973: #endif
974: for (j=0; j<n; j++) {
975: if (indexn[j] < 0) { idx++; continue;}
976: #if defined(PETSC_USE_DEBUG)
977: if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
978: #endif
979: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
980: }
981: }
982: } else {
983: for (i=0; i<m; i++) {
984: if (indexm[i] < 0) { idx += n; continue;}
985: #if defined(PETSC_USE_DEBUG)
986: if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1);
987: #endif
988: for (j=0; j<n; j++) {
989: if (indexn[j] < 0) { idx++; continue;}
990: #if defined(PETSC_USE_DEBUG)
991: if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1);
992: #endif
993: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
994: }
995: }
996: }
997: }
998: return(0);
999: }
1001: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1002: {
1003: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1004: PetscInt i,j;
1007: /* row-oriented output */
1008: for (i=0; i<m; i++) {
1009: if (indexm[i] < 0) {v += n;continue;}
1010: if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested larger than number rows %D",indexm[i],A->rmap->n);
1011: for (j=0; j<n; j++) {
1012: if (indexn[j] < 0) {v++; continue;}
1013: if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column %D requested larger than number columns %D",indexn[j],A->cmap->n);
1014: *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
1015: }
1016: }
1017: return(0);
1018: }
1020: /* -----------------------------------------------------------------*/
1022: static PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
1023: {
1024: Mat_SeqDense *a;
1026: PetscInt *scols,i,j,nz,header[4];
1027: int fd;
1028: PetscMPIInt size;
1029: PetscInt *rowlengths = 0,M,N,*cols,grows,gcols;
1030: PetscScalar *vals,*svals,*v,*w;
1031: MPI_Comm comm;
1034: /* force binary viewer to load .info file if it has not yet done so */
1035: PetscViewerSetUp(viewer);
1036: PetscObjectGetComm((PetscObject)viewer,&comm);
1037: MPI_Comm_size(comm,&size);
1038: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
1039: PetscViewerBinaryGetDescriptor(viewer,&fd);
1040: PetscBinaryRead(fd,header,4,PETSC_INT);
1041: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
1042: M = header[1]; N = header[2]; nz = header[3];
1044: /* set global size if not set already*/
1045: if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
1046: MatSetSizes(newmat,M,N,M,N);
1047: } else {
1048: /* if sizes and type are already set, check if the vector global sizes are correct */
1049: MatGetSize(newmat,&grows,&gcols);
1050: if (M != grows || N != gcols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,grows,gcols);
1051: }
1052: a = (Mat_SeqDense*)newmat->data;
1053: if (!a->user_alloc) {
1054: MatSeqDenseSetPreallocation(newmat,NULL);
1055: }
1057: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
1058: a = (Mat_SeqDense*)newmat->data;
1059: v = a->v;
1060: /* Allocate some temp space to read in the values and then flip them
1061: from row major to column major */
1062: PetscMalloc1(M*N > 0 ? M*N : 1,&w);
1063: /* read in nonzero values */
1064: PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
1065: /* now flip the values and store them in the matrix*/
1066: for (j=0; j<N; j++) {
1067: for (i=0; i<M; i++) {
1068: *v++ =w[i*N+j];
1069: }
1070: }
1071: PetscFree(w);
1072: } else {
1073: /* read row lengths */
1074: PetscMalloc1(M+1,&rowlengths);
1075: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
1077: a = (Mat_SeqDense*)newmat->data;
1078: v = a->v;
1080: /* read column indices and nonzeros */
1081: PetscMalloc1(nz+1,&scols);
1082: cols = scols;
1083: PetscBinaryRead(fd,cols,nz,PETSC_INT);
1084: PetscMalloc1(nz+1,&svals);
1085: vals = svals;
1086: PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);
1088: /* insert into matrix */
1089: for (i=0; i<M; i++) {
1090: for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
1091: svals += rowlengths[i]; scols += rowlengths[i];
1092: }
1093: PetscFree(vals);
1094: PetscFree(cols);
1095: PetscFree(rowlengths);
1096: }
1097: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1098: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1099: return(0);
1100: }
1102: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1103: {
1104: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1105: PetscErrorCode ierr;
1106: PetscInt i,j;
1107: const char *name;
1108: PetscScalar *v;
1109: PetscViewerFormat format;
1110: #if defined(PETSC_USE_COMPLEX)
1111: PetscBool allreal = PETSC_TRUE;
1112: #endif
1115: PetscViewerGetFormat(viewer,&format);
1116: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1117: return(0); /* do nothing for now */
1118: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1119: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1120: for (i=0; i<A->rmap->n; i++) {
1121: v = a->v + i;
1122: PetscViewerASCIIPrintf(viewer,"row %D:",i);
1123: for (j=0; j<A->cmap->n; j++) {
1124: #if defined(PETSC_USE_COMPLEX)
1125: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1126: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1127: } else if (PetscRealPart(*v)) {
1128: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1129: }
1130: #else
1131: if (*v) {
1132: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1133: }
1134: #endif
1135: v += a->lda;
1136: }
1137: PetscViewerASCIIPrintf(viewer,"\n");
1138: }
1139: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1140: } else {
1141: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1142: #if defined(PETSC_USE_COMPLEX)
1143: /* determine if matrix has all real values */
1144: v = a->v;
1145: for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1146: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1147: }
1148: #endif
1149: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1150: PetscObjectGetName((PetscObject)A,&name);
1151: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1152: PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1153: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1154: }
1156: for (i=0; i<A->rmap->n; i++) {
1157: v = a->v + i;
1158: for (j=0; j<A->cmap->n; j++) {
1159: #if defined(PETSC_USE_COMPLEX)
1160: if (allreal) {
1161: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1162: } else {
1163: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1164: }
1165: #else
1166: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1167: #endif
1168: v += a->lda;
1169: }
1170: PetscViewerASCIIPrintf(viewer,"\n");
1171: }
1172: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1173: PetscViewerASCIIPrintf(viewer,"];\n");
1174: }
1175: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1176: }
1177: PetscViewerFlush(viewer);
1178: return(0);
1179: }
1181: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1182: {
1183: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1184: PetscErrorCode ierr;
1185: int fd;
1186: PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1187: PetscScalar *v,*anonz,*vals;
1188: PetscViewerFormat format;
1191: PetscViewerBinaryGetDescriptor(viewer,&fd);
1193: PetscViewerGetFormat(viewer,&format);
1194: if (format == PETSC_VIEWER_NATIVE) {
1195: /* store the matrix as a dense matrix */
1196: PetscMalloc1(4,&col_lens);
1198: col_lens[0] = MAT_FILE_CLASSID;
1199: col_lens[1] = m;
1200: col_lens[2] = n;
1201: col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;
1203: PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1204: PetscFree(col_lens);
1206: /* write out matrix, by rows */
1207: PetscMalloc1(m*n+1,&vals);
1208: v = a->v;
1209: for (j=0; j<n; j++) {
1210: for (i=0; i<m; i++) {
1211: vals[j + i*n] = *v++;
1212: }
1213: }
1214: PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1215: PetscFree(vals);
1216: } else {
1217: PetscMalloc1(4+nz,&col_lens);
1219: col_lens[0] = MAT_FILE_CLASSID;
1220: col_lens[1] = m;
1221: col_lens[2] = n;
1222: col_lens[3] = nz;
1224: /* store lengths of each row and write (including header) to file */
1225: for (i=0; i<m; i++) col_lens[4+i] = n;
1226: PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);
1228: /* Possibly should write in smaller increments, not whole matrix at once? */
1229: /* store column indices (zero start index) */
1230: ict = 0;
1231: for (i=0; i<m; i++) {
1232: for (j=0; j<n; j++) col_lens[ict++] = j;
1233: }
1234: PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1235: PetscFree(col_lens);
1237: /* store nonzero values */
1238: PetscMalloc1(nz+1,&anonz);
1239: ict = 0;
1240: for (i=0; i<m; i++) {
1241: v = a->v + i;
1242: for (j=0; j<n; j++) {
1243: anonz[ict++] = *v; v += a->lda;
1244: }
1245: }
1246: PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1247: PetscFree(anonz);
1248: }
1249: return(0);
1250: }
1252: #include <petscdraw.h>
1253: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1254: {
1255: Mat A = (Mat) Aa;
1256: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1257: PetscErrorCode ierr;
1258: PetscInt m = A->rmap->n,n = A->cmap->n,i,j;
1259: int color = PETSC_DRAW_WHITE;
1260: PetscScalar *v = a->v;
1261: PetscViewer viewer;
1262: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1263: PetscViewerFormat format;
1266: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1267: PetscViewerGetFormat(viewer,&format);
1268: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
1270: /* Loop over matrix elements drawing boxes */
1272: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1273: PetscDrawCollectiveBegin(draw);
1274: /* Blue for negative and Red for positive */
1275: for (j = 0; j < n; j++) {
1276: x_l = j; x_r = x_l + 1.0;
1277: for (i = 0; i < m; i++) {
1278: y_l = m - i - 1.0;
1279: y_r = y_l + 1.0;
1280: if (PetscRealPart(v[j*m+i]) > 0.) {
1281: color = PETSC_DRAW_RED;
1282: } else if (PetscRealPart(v[j*m+i]) < 0.) {
1283: color = PETSC_DRAW_BLUE;
1284: } else {
1285: continue;
1286: }
1287: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1288: }
1289: }
1290: PetscDrawCollectiveEnd(draw);
1291: } else {
1292: /* use contour shading to indicate magnitude of values */
1293: /* first determine max of all nonzero values */
1294: PetscReal minv = 0.0, maxv = 0.0;
1295: PetscDraw popup;
1297: for (i=0; i < m*n; i++) {
1298: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1299: }
1300: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1301: PetscDrawGetPopup(draw,&popup);
1302: PetscDrawScalePopup(popup,minv,maxv);
1304: PetscDrawCollectiveBegin(draw);
1305: for (j=0; j<n; j++) {
1306: x_l = j;
1307: x_r = x_l + 1.0;
1308: for (i=0; i<m; i++) {
1309: y_l = m - i - 1.0;
1310: y_r = y_l + 1.0;
1311: color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1312: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1313: }
1314: }
1315: PetscDrawCollectiveEnd(draw);
1316: }
1317: return(0);
1318: }
1320: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1321: {
1322: PetscDraw draw;
1323: PetscBool isnull;
1324: PetscReal xr,yr,xl,yl,h,w;
1328: PetscViewerDrawGetDraw(viewer,0,&draw);
1329: PetscDrawIsNull(draw,&isnull);
1330: if (isnull) return(0);
1332: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1333: xr += w; yr += h; xl = -w; yl = -h;
1334: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1335: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1336: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1337: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1338: PetscDrawSave(draw);
1339: return(0);
1340: }
1342: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1343: {
1345: PetscBool iascii,isbinary,isdraw;
1348: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1349: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1350: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1352: if (iascii) {
1353: MatView_SeqDense_ASCII(A,viewer);
1354: } else if (isbinary) {
1355: MatView_SeqDense_Binary(A,viewer);
1356: } else if (isdraw) {
1357: MatView_SeqDense_Draw(A,viewer);
1358: }
1359: return(0);
1360: }
1362: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar array[])
1363: {
1364: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1367: a->unplacedarray = a->v;
1368: a->unplaced_user_alloc = a->user_alloc;
1369: a->v = (PetscScalar*) array;
1370: return(0);
1371: }
1373: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1374: {
1375: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1378: a->v = a->unplacedarray;
1379: a->user_alloc = a->unplaced_user_alloc;
1380: a->unplacedarray = NULL;
1381: return(0);
1382: }
1384: static PetscErrorCode MatDestroy_SeqDense(Mat mat)
1385: {
1386: Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
1390: #if defined(PETSC_USE_LOG)
1391: PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1392: #endif
1393: PetscFree(l->pivots);
1394: PetscFree(l->fwork);
1395: MatDestroy(&l->ptapwork);
1396: if (!l->user_alloc) {PetscFree(l->v);}
1397: PetscFree(mat->data);
1399: PetscObjectChangeTypeName((PetscObject)mat,0);
1400: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1401: PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1402: PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1403: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1404: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1405: #if defined(PETSC_HAVE_ELEMENTAL)
1406: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1407: #endif
1408: PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1409: PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1410: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1411: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1412: PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);
1413: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1414: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1415: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1416: return(0);
1417: }
1419: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1420: {
1421: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1423: PetscInt k,j,m,n,M;
1424: PetscScalar *v,tmp;
1427: v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1428: if (reuse == MAT_INPLACE_MATRIX) { /* in place transpose */
1429: if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place");
1430: else {
1431: for (j=0; j<m; j++) {
1432: for (k=0; k<j; k++) {
1433: tmp = v[j + k*M];
1434: v[j + k*M] = v[k + j*M];
1435: v[k + j*M] = tmp;
1436: }
1437: }
1438: }
1439: } else { /* out-of-place transpose */
1440: Mat tmat;
1441: Mat_SeqDense *tmatd;
1442: PetscScalar *v2;
1443: PetscInt M2;
1445: if (reuse == MAT_INITIAL_MATRIX) {
1446: MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1447: MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1448: MatSetType(tmat,((PetscObject)A)->type_name);
1449: MatSeqDenseSetPreallocation(tmat,NULL);
1450: } else {
1451: tmat = *matout;
1452: }
1453: tmatd = (Mat_SeqDense*)tmat->data;
1454: v = mat->v; v2 = tmatd->v; M2 = tmatd->lda;
1455: for (j=0; j<n; j++) {
1456: for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1457: }
1458: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1459: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1460: *matout = tmat;
1461: }
1462: return(0);
1463: }
1465: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg)
1466: {
1467: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1468: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1469: PetscInt i,j;
1470: PetscScalar *v1,*v2;
1473: if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1474: if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1475: for (i=0; i<A1->rmap->n; i++) {
1476: v1 = mat1->v+i; v2 = mat2->v+i;
1477: for (j=0; j<A1->cmap->n; j++) {
1478: if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1479: v1 += mat1->lda; v2 += mat2->lda;
1480: }
1481: }
1482: *flg = PETSC_TRUE;
1483: return(0);
1484: }
1486: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1487: {
1488: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1490: PetscInt i,n,len;
1491: PetscScalar *x,zero = 0.0;
1494: VecSet(v,zero);
1495: VecGetSize(v,&n);
1496: VecGetArray(v,&x);
1497: len = PetscMin(A->rmap->n,A->cmap->n);
1498: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1499: for (i=0; i<len; i++) {
1500: x[i] = mat->v[i*mat->lda + i];
1501: }
1502: VecRestoreArray(v,&x);
1503: return(0);
1504: }
1506: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1507: {
1508: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1509: const PetscScalar *l,*r;
1510: PetscScalar x,*v;
1511: PetscErrorCode ierr;
1512: PetscInt i,j,m = A->rmap->n,n = A->cmap->n;
1515: if (ll) {
1516: VecGetSize(ll,&m);
1517: VecGetArrayRead(ll,&l);
1518: if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1519: for (i=0; i<m; i++) {
1520: x = l[i];
1521: v = mat->v + i;
1522: for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1523: }
1524: VecRestoreArrayRead(ll,&l);
1525: PetscLogFlops(1.0*n*m);
1526: }
1527: if (rr) {
1528: VecGetSize(rr,&n);
1529: VecGetArrayRead(rr,&r);
1530: if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1531: for (i=0; i<n; i++) {
1532: x = r[i];
1533: v = mat->v + i*mat->lda;
1534: for (j=0; j<m; j++) (*v++) *= x;
1535: }
1536: VecRestoreArrayRead(rr,&r);
1537: PetscLogFlops(1.0*n*m);
1538: }
1539: return(0);
1540: }
1542: static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1543: {
1544: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1545: PetscScalar *v = mat->v;
1546: PetscReal sum = 0.0;
1547: PetscInt lda =mat->lda,m=A->rmap->n,i,j;
1551: if (type == NORM_FROBENIUS) {
1552: if (lda>m) {
1553: for (j=0; j<A->cmap->n; j++) {
1554: v = mat->v+j*lda;
1555: for (i=0; i<m; i++) {
1556: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1557: }
1558: }
1559: } else {
1560: #if defined(PETSC_USE_REAL___FP16)
1561: PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1562: *nrm = BLASnrm2_(&cnt,v,&one);
1563: }
1564: #else
1565: for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1566: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1567: }
1568: }
1569: *nrm = PetscSqrtReal(sum);
1570: #endif
1571: PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1572: } else if (type == NORM_1) {
1573: *nrm = 0.0;
1574: for (j=0; j<A->cmap->n; j++) {
1575: v = mat->v + j*mat->lda;
1576: sum = 0.0;
1577: for (i=0; i<A->rmap->n; i++) {
1578: sum += PetscAbsScalar(*v); v++;
1579: }
1580: if (sum > *nrm) *nrm = sum;
1581: }
1582: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1583: } else if (type == NORM_INFINITY) {
1584: *nrm = 0.0;
1585: for (j=0; j<A->rmap->n; j++) {
1586: v = mat->v + j;
1587: sum = 0.0;
1588: for (i=0; i<A->cmap->n; i++) {
1589: sum += PetscAbsScalar(*v); v += mat->lda;
1590: }
1591: if (sum > *nrm) *nrm = sum;
1592: }
1593: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1594: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1595: return(0);
1596: }
1598: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1599: {
1600: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
1604: switch (op) {
1605: case MAT_ROW_ORIENTED:
1606: aij->roworiented = flg;
1607: break;
1608: case MAT_NEW_NONZERO_LOCATIONS:
1609: case MAT_NEW_NONZERO_LOCATION_ERR:
1610: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1611: case MAT_NEW_DIAGONALS:
1612: case MAT_KEEP_NONZERO_PATTERN:
1613: case MAT_IGNORE_OFF_PROC_ENTRIES:
1614: case MAT_USE_HASH_TABLE:
1615: case MAT_IGNORE_ZERO_ENTRIES:
1616: case MAT_IGNORE_LOWER_TRIANGULAR:
1617: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1618: break;
1619: case MAT_SPD:
1620: case MAT_SYMMETRIC:
1621: case MAT_STRUCTURALLY_SYMMETRIC:
1622: case MAT_HERMITIAN:
1623: case MAT_SYMMETRY_ETERNAL:
1624: /* These options are handled directly by MatSetOption() */
1625: break;
1626: default:
1627: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1628: }
1629: return(0);
1630: }
1632: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1633: {
1634: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1636: PetscInt lda=l->lda,m=A->rmap->n,j;
1639: if (lda>m) {
1640: for (j=0; j<A->cmap->n; j++) {
1641: PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1642: }
1643: } else {
1644: PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1645: }
1646: return(0);
1647: }
1649: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1650: {
1651: PetscErrorCode ierr;
1652: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1653: PetscInt m = l->lda, n = A->cmap->n, i,j;
1654: PetscScalar *slot,*bb;
1655: const PetscScalar *xx;
1658: #if defined(PETSC_USE_DEBUG)
1659: for (i=0; i<N; i++) {
1660: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1661: if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n);
1662: }
1663: #endif
1665: /* fix right hand side if needed */
1666: if (x && b) {
1667: VecGetArrayRead(x,&xx);
1668: VecGetArray(b,&bb);
1669: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1670: VecRestoreArrayRead(x,&xx);
1671: VecRestoreArray(b,&bb);
1672: }
1674: for (i=0; i<N; i++) {
1675: slot = l->v + rows[i];
1676: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1677: }
1678: if (diag != 0.0) {
1679: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1680: for (i=0; i<N; i++) {
1681: slot = l->v + (m+1)*rows[i];
1682: *slot = diag;
1683: }
1684: }
1685: return(0);
1686: }
1688: static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1689: {
1690: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1693: if (mat->lda != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get array for Dense matrices with LDA different from number of rows");
1694: *array = mat->v;
1695: return(0);
1696: }
1698: static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1699: {
1701: *array = 0; /* user cannot accidently use the array later */
1702: return(0);
1703: }
1705: /*@C
1706: MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored
1708: Not Collective
1710: Input Parameter:
1711: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1713: Output Parameter:
1714: . array - pointer to the data
1716: Level: intermediate
1718: .seealso: MatDenseRestoreArray()
1719: @*/
1720: PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array)
1721: {
1725: PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1726: return(0);
1727: }
1729: /*@C
1730: MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
1732: Not Collective
1734: Input Parameters:
1735: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1736: . array - pointer to the data
1738: Level: intermediate
1740: .seealso: MatDenseGetArray()
1741: @*/
1742: PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array)
1743: {
1747: PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1748: return(0);
1749: }
1751: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1752: {
1753: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1755: PetscInt i,j,nrows,ncols;
1756: const PetscInt *irow,*icol;
1757: PetscScalar *av,*bv,*v = mat->v;
1758: Mat newmat;
1761: ISGetIndices(isrow,&irow);
1762: ISGetIndices(iscol,&icol);
1763: ISGetLocalSize(isrow,&nrows);
1764: ISGetLocalSize(iscol,&ncols);
1766: /* Check submatrixcall */
1767: if (scall == MAT_REUSE_MATRIX) {
1768: PetscInt n_cols,n_rows;
1769: MatGetSize(*B,&n_rows,&n_cols);
1770: if (n_rows != nrows || n_cols != ncols) {
1771: /* resize the result matrix to match number of requested rows/columns */
1772: MatSetSizes(*B,nrows,ncols,nrows,ncols);
1773: }
1774: newmat = *B;
1775: } else {
1776: /* Create and fill new matrix */
1777: MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1778: MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1779: MatSetType(newmat,((PetscObject)A)->type_name);
1780: MatSeqDenseSetPreallocation(newmat,NULL);
1781: }
1783: /* Now extract the data pointers and do the copy,column at a time */
1784: bv = ((Mat_SeqDense*)newmat->data)->v;
1786: for (i=0; i<ncols; i++) {
1787: av = v + mat->lda*icol[i];
1788: for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1789: }
1791: /* Assemble the matrices so that the correct flags are set */
1792: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1793: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1795: /* Free work space */
1796: ISRestoreIndices(isrow,&irow);
1797: ISRestoreIndices(iscol,&icol);
1798: *B = newmat;
1799: return(0);
1800: }
1802: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
1803: {
1805: PetscInt i;
1808: if (scall == MAT_INITIAL_MATRIX) {
1809: PetscCalloc1(n+1,B);
1810: }
1812: for (i=0; i<n; i++) {
1813: MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1814: }
1815: return(0);
1816: }
1818: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1819: {
1821: return(0);
1822: }
1824: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1825: {
1827: return(0);
1828: }
1830: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1831: {
1832: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1834: PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;
1837: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1838: if (A->ops->copy != B->ops->copy) {
1839: MatCopy_Basic(A,B,str);
1840: return(0);
1841: }
1842: if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1843: if (lda1>m || lda2>m) {
1844: for (j=0; j<n; j++) {
1845: PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1846: }
1847: } else {
1848: PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1849: }
1850: PetscObjectStateIncrease((PetscObject)B);
1851: return(0);
1852: }
1854: static PetscErrorCode MatSetUp_SeqDense(Mat A)
1855: {
1859: MatSeqDenseSetPreallocation(A,0);
1860: return(0);
1861: }
1863: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1864: {
1865: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1866: PetscInt i,nz = A->rmap->n*A->cmap->n;
1867: PetscScalar *aa = a->v;
1870: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1871: return(0);
1872: }
1874: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1875: {
1876: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1877: PetscInt i,nz = A->rmap->n*A->cmap->n;
1878: PetscScalar *aa = a->v;
1881: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1882: return(0);
1883: }
1885: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1886: {
1887: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1888: PetscInt i,nz = A->rmap->n*A->cmap->n;
1889: PetscScalar *aa = a->v;
1892: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1893: return(0);
1894: }
1896: /* ----------------------------------------------------------------*/
1897: PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1898: {
1902: if (scall == MAT_INITIAL_MATRIX) {
1903: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1904: MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1905: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1906: }
1907: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1908: MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
1909: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1910: return(0);
1911: }
1913: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1914: {
1916: PetscInt m=A->rmap->n,n=B->cmap->n;
1917: Mat Cmat;
1920: 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);
1921: MatCreate(PETSC_COMM_SELF,&Cmat);
1922: MatSetSizes(Cmat,m,n,m,n);
1923: MatSetType(Cmat,MATSEQDENSE);
1924: MatSeqDenseSetPreallocation(Cmat,NULL);
1926: *C = Cmat;
1927: return(0);
1928: }
1930: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1931: {
1932: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1933: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
1934: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
1935: PetscBLASInt m,n,k;
1936: PetscScalar _DOne=1.0,_DZero=0.0;
1938: PetscBool flg;
1941: PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
1942: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Second matrix must be dense");
1944: /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/
1945: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
1946: if (flg) {
1947: C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
1948: (*C->ops->matmultnumeric)(A,B,C);
1949: return(0);
1950: }
1952: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
1953: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense");
1954: PetscBLASIntCast(C->rmap->n,&m);
1955: PetscBLASIntCast(C->cmap->n,&n);
1956: PetscBLASIntCast(A->cmap->n,&k);
1957: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
1958: return(0);
1959: }
1961: PetscErrorCode MatMatTransposeMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1962: {
1966: if (scall == MAT_INITIAL_MATRIX) {
1967: PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);
1968: MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
1969: PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);
1970: }
1971: PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);
1972: MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);
1973: PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);
1974: return(0);
1975: }
1977: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1978: {
1980: PetscInt m=A->rmap->n,n=B->rmap->n;
1981: Mat Cmat;
1984: if (A->cmap->n != B->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->cmap->n %d\n",A->cmap->n,B->cmap->n);
1985: MatCreate(PETSC_COMM_SELF,&Cmat);
1986: MatSetSizes(Cmat,m,n,m,n);
1987: MatSetType(Cmat,MATSEQDENSE);
1988: MatSeqDenseSetPreallocation(Cmat,NULL);
1990: Cmat->assembled = PETSC_TRUE;
1992: *C = Cmat;
1993: return(0);
1994: }
1996: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
1997: {
1998: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1999: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2000: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2001: PetscBLASInt m,n,k;
2002: PetscScalar _DOne=1.0,_DZero=0.0;
2006: PetscBLASIntCast(A->rmap->n,&m);
2007: PetscBLASIntCast(B->rmap->n,&n);
2008: PetscBLASIntCast(A->cmap->n,&k);
2009: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2010: return(0);
2011: }
2013: PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2014: {
2018: if (scall == MAT_INITIAL_MATRIX) {
2019: PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
2020: MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2021: PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
2022: }
2023: PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
2024: MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2025: PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
2026: return(0);
2027: }
2029: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2030: {
2032: PetscInt m=A->cmap->n,n=B->cmap->n;
2033: Mat Cmat;
2036: if (A->rmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->rmap->n %d != B->rmap->n %d\n",A->rmap->n,B->rmap->n);
2037: MatCreate(PETSC_COMM_SELF,&Cmat);
2038: MatSetSizes(Cmat,m,n,m,n);
2039: MatSetType(Cmat,MATSEQDENSE);
2040: MatSeqDenseSetPreallocation(Cmat,NULL);
2042: Cmat->assembled = PETSC_TRUE;
2044: *C = Cmat;
2045: return(0);
2046: }
2048: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2049: {
2050: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2051: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2052: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2053: PetscBLASInt m,n,k;
2054: PetscScalar _DOne=1.0,_DZero=0.0;
2058: PetscBLASIntCast(C->rmap->n,&m);
2059: PetscBLASIntCast(C->cmap->n,&n);
2060: PetscBLASIntCast(A->rmap->n,&k);
2061: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2062: return(0);
2063: }
2065: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2066: {
2067: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2069: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2070: PetscScalar *x;
2071: MatScalar *aa = a->v;
2074: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2076: VecSet(v,0.0);
2077: VecGetArray(v,&x);
2078: VecGetLocalSize(v,&p);
2079: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2080: for (i=0; i<m; i++) {
2081: x[i] = aa[i]; if (idx) idx[i] = 0;
2082: for (j=1; j<n; j++) {
2083: if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2084: }
2085: }
2086: VecRestoreArray(v,&x);
2087: return(0);
2088: }
2090: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2091: {
2092: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2094: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2095: PetscScalar *x;
2096: PetscReal atmp;
2097: MatScalar *aa = a->v;
2100: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2102: VecSet(v,0.0);
2103: VecGetArray(v,&x);
2104: VecGetLocalSize(v,&p);
2105: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2106: for (i=0; i<m; i++) {
2107: x[i] = PetscAbsScalar(aa[i]);
2108: for (j=1; j<n; j++) {
2109: atmp = PetscAbsScalar(aa[i+m*j]);
2110: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2111: }
2112: }
2113: VecRestoreArray(v,&x);
2114: return(0);
2115: }
2117: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2118: {
2119: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2121: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2122: PetscScalar *x;
2123: MatScalar *aa = a->v;
2126: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2128: VecSet(v,0.0);
2129: VecGetArray(v,&x);
2130: VecGetLocalSize(v,&p);
2131: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2132: for (i=0; i<m; i++) {
2133: x[i] = aa[i]; if (idx) idx[i] = 0;
2134: for (j=1; j<n; j++) {
2135: if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2136: }
2137: }
2138: VecRestoreArray(v,&x);
2139: return(0);
2140: }
2142: static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2143: {
2144: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2146: PetscScalar *x;
2149: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2151: VecGetArray(v,&x);
2152: PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
2153: VecRestoreArray(v,&x);
2154: return(0);
2155: }
2158: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2159: {
2161: PetscInt i,j,m,n;
2162: PetscScalar *a;
2165: MatGetSize(A,&m,&n);
2166: PetscMemzero(norms,n*sizeof(PetscReal));
2167: MatDenseGetArray(A,&a);
2168: if (type == NORM_2) {
2169: for (i=0; i<n; i++) {
2170: for (j=0; j<m; j++) {
2171: norms[i] += PetscAbsScalar(a[j]*a[j]);
2172: }
2173: a += m;
2174: }
2175: } else if (type == NORM_1) {
2176: for (i=0; i<n; i++) {
2177: for (j=0; j<m; j++) {
2178: norms[i] += PetscAbsScalar(a[j]);
2179: }
2180: a += m;
2181: }
2182: } else if (type == NORM_INFINITY) {
2183: for (i=0; i<n; i++) {
2184: for (j=0; j<m; j++) {
2185: norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2186: }
2187: a += m;
2188: }
2189: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2190: MatDenseRestoreArray(A,&a);
2191: if (type == NORM_2) {
2192: for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2193: }
2194: return(0);
2195: }
2197: static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2198: {
2200: PetscScalar *a;
2201: PetscInt m,n,i;
2204: MatGetSize(x,&m,&n);
2205: MatDenseGetArray(x,&a);
2206: for (i=0; i<m*n; i++) {
2207: PetscRandomGetValue(rctx,a+i);
2208: }
2209: MatDenseRestoreArray(x,&a);
2210: return(0);
2211: }
2213: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d)
2214: {
2216: *missing = PETSC_FALSE;
2217: return(0);
2218: }
2221: /* -------------------------------------------------------------------*/
2222: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2223: MatGetRow_SeqDense,
2224: MatRestoreRow_SeqDense,
2225: MatMult_SeqDense,
2226: /* 4*/ MatMultAdd_SeqDense,
2227: MatMultTranspose_SeqDense,
2228: MatMultTransposeAdd_SeqDense,
2229: 0,
2230: 0,
2231: 0,
2232: /* 10*/ 0,
2233: MatLUFactor_SeqDense,
2234: MatCholeskyFactor_SeqDense,
2235: MatSOR_SeqDense,
2236: MatTranspose_SeqDense,
2237: /* 15*/ MatGetInfo_SeqDense,
2238: MatEqual_SeqDense,
2239: MatGetDiagonal_SeqDense,
2240: MatDiagonalScale_SeqDense,
2241: MatNorm_SeqDense,
2242: /* 20*/ MatAssemblyBegin_SeqDense,
2243: MatAssemblyEnd_SeqDense,
2244: MatSetOption_SeqDense,
2245: MatZeroEntries_SeqDense,
2246: /* 24*/ MatZeroRows_SeqDense,
2247: 0,
2248: 0,
2249: 0,
2250: 0,
2251: /* 29*/ MatSetUp_SeqDense,
2252: 0,
2253: 0,
2254: 0,
2255: 0,
2256: /* 34*/ MatDuplicate_SeqDense,
2257: 0,
2258: 0,
2259: 0,
2260: 0,
2261: /* 39*/ MatAXPY_SeqDense,
2262: MatCreateSubMatrices_SeqDense,
2263: 0,
2264: MatGetValues_SeqDense,
2265: MatCopy_SeqDense,
2266: /* 44*/ MatGetRowMax_SeqDense,
2267: MatScale_SeqDense,
2268: MatShift_Basic,
2269: 0,
2270: MatZeroRowsColumns_SeqDense,
2271: /* 49*/ MatSetRandom_SeqDense,
2272: 0,
2273: 0,
2274: 0,
2275: 0,
2276: /* 54*/ 0,
2277: 0,
2278: 0,
2279: 0,
2280: 0,
2281: /* 59*/ 0,
2282: MatDestroy_SeqDense,
2283: MatView_SeqDense,
2284: 0,
2285: 0,
2286: /* 64*/ 0,
2287: 0,
2288: 0,
2289: 0,
2290: 0,
2291: /* 69*/ MatGetRowMaxAbs_SeqDense,
2292: 0,
2293: 0,
2294: 0,
2295: 0,
2296: /* 74*/ 0,
2297: 0,
2298: 0,
2299: 0,
2300: 0,
2301: /* 79*/ 0,
2302: 0,
2303: 0,
2304: 0,
2305: /* 83*/ MatLoad_SeqDense,
2306: 0,
2307: MatIsHermitian_SeqDense,
2308: 0,
2309: 0,
2310: 0,
2311: /* 89*/ MatMatMult_SeqDense_SeqDense,
2312: MatMatMultSymbolic_SeqDense_SeqDense,
2313: MatMatMultNumeric_SeqDense_SeqDense,
2314: MatPtAP_SeqDense_SeqDense,
2315: MatPtAPSymbolic_SeqDense_SeqDense,
2316: /* 94*/ MatPtAPNumeric_SeqDense_SeqDense,
2317: MatMatTransposeMult_SeqDense_SeqDense,
2318: MatMatTransposeMultSymbolic_SeqDense_SeqDense,
2319: MatMatTransposeMultNumeric_SeqDense_SeqDense,
2320: 0,
2321: /* 99*/ 0,
2322: 0,
2323: 0,
2324: MatConjugate_SeqDense,
2325: 0,
2326: /*104*/ 0,
2327: MatRealPart_SeqDense,
2328: MatImaginaryPart_SeqDense,
2329: 0,
2330: 0,
2331: /*109*/ 0,
2332: 0,
2333: MatGetRowMin_SeqDense,
2334: MatGetColumnVector_SeqDense,
2335: MatMissingDiagonal_SeqDense,
2336: /*114*/ 0,
2337: 0,
2338: 0,
2339: 0,
2340: 0,
2341: /*119*/ 0,
2342: 0,
2343: 0,
2344: 0,
2345: 0,
2346: /*124*/ 0,
2347: MatGetColumnNorms_SeqDense,
2348: 0,
2349: 0,
2350: 0,
2351: /*129*/ 0,
2352: MatTransposeMatMult_SeqDense_SeqDense,
2353: MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2354: MatTransposeMatMultNumeric_SeqDense_SeqDense,
2355: 0,
2356: /*134*/ 0,
2357: 0,
2358: 0,
2359: 0,
2360: 0,
2361: /*139*/ 0,
2362: 0,
2363: 0
2364: };
2366: /*@C
2367: MatCreateSeqDense - Creates a sequential dense matrix that
2368: is stored in column major order (the usual Fortran 77 manner). Many
2369: of the matrix operations use the BLAS and LAPACK routines.
2371: Collective on MPI_Comm
2373: Input Parameters:
2374: + comm - MPI communicator, set to PETSC_COMM_SELF
2375: . m - number of rows
2376: . n - number of columns
2377: - data - optional location of matrix data in column major order. Set data=NULL for PETSc
2378: to control all matrix memory allocation.
2380: Output Parameter:
2381: . A - the matrix
2383: Notes:
2384: The data input variable is intended primarily for Fortran programmers
2385: who wish to allocate their own matrix memory space. Most users should
2386: set data=NULL.
2388: Level: intermediate
2390: .keywords: dense, matrix, LAPACK, BLAS
2392: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2393: @*/
2394: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2395: {
2399: MatCreate(comm,A);
2400: MatSetSizes(*A,m,n,m,n);
2401: MatSetType(*A,MATSEQDENSE);
2402: MatSeqDenseSetPreallocation(*A,data);
2403: return(0);
2404: }
2406: /*@C
2407: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
2409: Collective on MPI_Comm
2411: Input Parameters:
2412: + B - the matrix
2413: - data - the array (or NULL)
2415: Notes:
2416: The data input variable is intended primarily for Fortran programmers
2417: who wish to allocate their own matrix memory space. Most users should
2418: need not call this routine.
2420: Level: intermediate
2422: .keywords: dense, matrix, LAPACK, BLAS
2424: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA()
2426: @*/
2427: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2428: {
2432: PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2433: return(0);
2434: }
2436: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2437: {
2438: Mat_SeqDense *b;
2442: B->preallocated = PETSC_TRUE;
2444: PetscLayoutSetUp(B->rmap);
2445: PetscLayoutSetUp(B->cmap);
2447: b = (Mat_SeqDense*)B->data;
2448: b->Mmax = B->rmap->n;
2449: b->Nmax = B->cmap->n;
2450: if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;
2452: PetscIntMultError(b->lda,b->Nmax,NULL);
2453: if (!data) { /* petsc-allocated storage */
2454: if (!b->user_alloc) { PetscFree(b->v); }
2455: PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2456: PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));
2458: b->user_alloc = PETSC_FALSE;
2459: } else { /* user-allocated storage */
2460: if (!b->user_alloc) { PetscFree(b->v); }
2461: b->v = data;
2462: b->user_alloc = PETSC_TRUE;
2463: }
2464: B->assembled = PETSC_TRUE;
2465: return(0);
2466: }
2468: #if defined(PETSC_HAVE_ELEMENTAL)
2469: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2470: {
2471: Mat mat_elemental;
2473: PetscScalar *array,*v_colwise;
2474: PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;
2477: PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2478: MatDenseGetArray(A,&array);
2479: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2480: k = 0;
2481: for (j=0; j<N; j++) {
2482: cols[j] = j;
2483: for (i=0; i<M; i++) {
2484: v_colwise[j*M+i] = array[k++];
2485: }
2486: }
2487: for (i=0; i<M; i++) {
2488: rows[i] = i;
2489: }
2490: MatDenseRestoreArray(A,&array);
2492: MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2493: MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2494: MatSetType(mat_elemental,MATELEMENTAL);
2495: MatSetUp(mat_elemental);
2497: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2498: MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2499: MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2500: MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2501: PetscFree3(v_colwise,rows,cols);
2503: if (reuse == MAT_INPLACE_MATRIX) {
2504: MatHeaderReplace(A,&mat_elemental);
2505: } else {
2506: *newmat = mat_elemental;
2507: }
2508: return(0);
2509: }
2510: #endif
2512: /*@C
2513: MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array
2515: Input parameter:
2516: + A - the matrix
2517: - lda - the leading dimension
2519: Notes:
2520: This routine is to be used in conjunction with MatSeqDenseSetPreallocation();
2521: it asserts that the preallocation has a leading dimension (the LDA parameter
2522: of Blas and Lapack fame) larger than M, the first dimension of the matrix.
2524: Level: intermediate
2526: .keywords: dense, matrix, LAPACK, BLAS
2528: .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize()
2530: @*/
2531: PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda)
2532: {
2533: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2536: if (lda < B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap->n);
2537: b->lda = lda;
2538: b->changelda = PETSC_FALSE;
2539: b->Mmax = PetscMax(b->Mmax,lda);
2540: return(0);
2541: }
2543: /*MC
2544: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
2546: Options Database Keys:
2547: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
2549: Level: beginner
2551: .seealso: MatCreateSeqDense()
2553: M*/
2555: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2556: {
2557: Mat_SeqDense *b;
2559: PetscMPIInt size;
2562: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
2563: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
2565: PetscNewLog(B,&b);
2566: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2567: B->data = (void*)b;
2569: b->roworiented = PETSC_TRUE;
2571: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2572: PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
2573: PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
2574: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2575: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2576: #if defined(PETSC_HAVE_ELEMENTAL)
2577: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2578: #endif
2579: PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2580: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2581: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2582: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2583: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);
2584: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2585: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2586: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2587: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);
2588: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2589: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2590: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2591: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);
2593: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2594: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2595: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2596: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2597: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2598: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2600: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2601: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2602: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2603: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2604: return(0);
2605: }