Actual source code: dense.c
petsc-3.11.4 2019-09-28
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: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
55: PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
56: PetscFPTrapPop();
57: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0); /* TODO CHECK FLOPS */
58: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
59: if (A->spd) {
60: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
61: PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info));
62: PetscFPTrapPop();
63: MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
64: #if defined(PETSC_USE_COMPLEX)
65: } else if (A->hermitian) {
66: if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
67: if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present");
68: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
69: PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
70: PetscFPTrapPop();
71: MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);
72: #endif
73: } else { /* symmetric case */
74: if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present");
75: if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present");
76: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
77: PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
78: PetscFPTrapPop();
79: MatSeqDenseSymmetrize_Private(A,PETSC_FALSE);
80: }
81: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad Inversion: zero pivot in row %D",(PetscInt)info-1);
82: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0); /* TODO CHECK FLOPS */
83: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
84: #endif
86: A->ops->solve = NULL;
87: A->ops->matsolve = NULL;
88: A->ops->solvetranspose = NULL;
89: A->ops->matsolvetranspose = NULL;
90: A->ops->solveadd = NULL;
91: A->ops->solvetransposeadd = NULL;
92: A->factortype = MAT_FACTOR_NONE;
93: PetscFree(A->solvertype);
94: return(0);
95: }
97: PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
98: {
99: PetscErrorCode ierr;
100: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
101: PetscInt m = l->lda, n = A->cmap->n,r = A->rmap->n, i,j;
102: PetscScalar *slot,*bb;
103: const PetscScalar *xx;
106: #if defined(PETSC_USE_DEBUG)
107: for (i=0; i<N; i++) {
108: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
109: 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);
110: 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);
111: }
112: #endif
114: /* fix right hand side if needed */
115: if (x && b) {
116: Vec xt;
118: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
119: VecDuplicate(x,&xt);
120: VecCopy(x,xt);
121: VecScale(xt,-1.0);
122: MatMultAdd(A,xt,b,b);
123: VecDestroy(&xt);
124: VecGetArrayRead(x,&xx);
125: VecGetArray(b,&bb);
126: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
127: VecRestoreArrayRead(x,&xx);
128: VecRestoreArray(b,&bb);
129: }
131: for (i=0; i<N; i++) {
132: slot = l->v + rows[i]*m;
133: PetscMemzero(slot,r*sizeof(PetscScalar));
134: }
135: for (i=0; i<N; i++) {
136: slot = l->v + rows[i];
137: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
138: }
139: if (diag != 0.0) {
140: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
141: for (i=0; i<N; i++) {
142: slot = l->v + (m+1)*rows[i];
143: *slot = diag;
144: }
145: }
146: return(0);
147: }
149: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
150: {
151: Mat_SeqDense *c = (Mat_SeqDense*)(C->data);
155: MatMatMultNumeric_SeqDense_SeqDense(A,P,c->ptapwork);
156: MatTransposeMatMultNumeric_SeqDense_SeqDense(P,c->ptapwork,C);
157: return(0);
158: }
160: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat *C)
161: {
162: Mat_SeqDense *c;
166: MatCreateSeqDense(PetscObjectComm((PetscObject)A),P->cmap->N,P->cmap->N,NULL,C);
167: c = (Mat_SeqDense*)((*C)->data);
168: MatCreateSeqDense(PetscObjectComm((PetscObject)A),A->rmap->N,P->cmap->N,NULL,&c->ptapwork);
169: return(0);
170: }
172: PETSC_INTERN PetscErrorCode MatPtAP_SeqDense_SeqDense(Mat A,Mat P,MatReuse reuse,PetscReal fill,Mat *C)
173: {
177: if (reuse == MAT_INITIAL_MATRIX) {
178: MatPtAPSymbolic_SeqDense_SeqDense(A,P,fill,C);
179: }
180: PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
181: (*(*C)->ops->ptapnumeric)(A,P,*C);
182: PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
183: return(0);
184: }
186: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
187: {
188: Mat B = NULL;
189: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
190: Mat_SeqDense *b;
192: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
193: MatScalar *av=a->a;
194: PetscBool isseqdense;
197: if (reuse == MAT_REUSE_MATRIX) {
198: PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);
199: if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type_name);
200: }
201: if (reuse != MAT_REUSE_MATRIX) {
202: MatCreate(PetscObjectComm((PetscObject)A),&B);
203: MatSetSizes(B,m,n,m,n);
204: MatSetType(B,MATSEQDENSE);
205: MatSeqDenseSetPreallocation(B,NULL);
206: b = (Mat_SeqDense*)(B->data);
207: } else {
208: b = (Mat_SeqDense*)((*newmat)->data);
209: PetscMemzero(b->v,m*n*sizeof(PetscScalar));
210: }
211: for (i=0; i<m; i++) {
212: PetscInt j;
213: for (j=0;j<ai[1]-ai[0];j++) {
214: b->v[*aj*m+i] = *av;
215: aj++;
216: av++;
217: }
218: ai++;
219: }
221: if (reuse == MAT_INPLACE_MATRIX) {
222: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
223: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
224: MatHeaderReplace(A,&B);
225: } else {
226: if (B) *newmat = B;
227: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
228: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
229: }
230: return(0);
231: }
233: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
234: {
235: Mat B;
236: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
238: PetscInt i, j;
239: PetscInt *rows, *nnz;
240: MatScalar *aa = a->v, *vals;
243: MatCreate(PetscObjectComm((PetscObject)A),&B);
244: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
245: MatSetType(B,MATSEQAIJ);
246: PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
247: for (j=0; j<A->cmap->n; j++) {
248: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
249: aa += a->lda;
250: }
251: MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
252: aa = a->v;
253: for (j=0; j<A->cmap->n; j++) {
254: PetscInt numRows = 0;
255: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
256: MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
257: aa += a->lda;
258: }
259: PetscFree3(rows,nnz,vals);
260: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
261: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
263: if (reuse == MAT_INPLACE_MATRIX) {
264: MatHeaderReplace(A,&B);
265: } else {
266: *newmat = B;
267: }
268: return(0);
269: }
271: static PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
272: {
273: Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
274: PetscScalar oalpha = alpha;
275: PetscInt j;
276: PetscBLASInt N,m,ldax,lday,one = 1;
280: PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
281: PetscBLASIntCast(X->rmap->n,&m);
282: PetscBLASIntCast(x->lda,&ldax);
283: PetscBLASIntCast(y->lda,&lday);
284: if (ldax>m || lday>m) {
285: for (j=0; j<X->cmap->n; j++) {
286: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one));
287: }
288: } else {
289: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one));
290: }
291: PetscObjectStateIncrease((PetscObject)Y);
292: PetscLogFlops(PetscMax(2*N-1,0));
293: return(0);
294: }
296: static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
297: {
298: PetscInt N = A->rmap->n*A->cmap->n;
301: info->block_size = 1.0;
302: info->nz_allocated = (double)N;
303: info->nz_used = (double)N;
304: info->nz_unneeded = (double)0;
305: info->assemblies = (double)A->num_ass;
306: info->mallocs = 0;
307: info->memory = ((PetscObject)A)->mem;
308: info->fill_ratio_given = 0;
309: info->fill_ratio_needed = 0;
310: info->factor_mallocs = 0;
311: return(0);
312: }
314: static PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
315: {
316: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
317: PetscScalar oalpha = alpha;
319: PetscBLASInt one = 1,j,nz,lda;
322: PetscBLASIntCast(a->lda,&lda);
323: if (lda>A->rmap->n) {
324: PetscBLASIntCast(A->rmap->n,&nz);
325: for (j=0; j<A->cmap->n; j++) {
326: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
327: }
328: } else {
329: PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
330: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
331: }
332: PetscLogFlops(nz);
333: return(0);
334: }
336: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl)
337: {
338: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
339: PetscInt i,j,m = A->rmap->n,N;
340: PetscScalar *v = a->v;
343: *fl = PETSC_FALSE;
344: if (A->rmap->n != A->cmap->n) return(0);
345: N = a->lda;
347: for (i=0; i<m; i++) {
348: for (j=i+1; j<m; j++) {
349: if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
350: }
351: }
352: *fl = PETSC_TRUE;
353: return(0);
354: }
356: static PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
357: {
358: Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l;
360: PetscInt lda = (PetscInt)mat->lda,j,m;
363: PetscLayoutReference(A->rmap,&newi->rmap);
364: PetscLayoutReference(A->cmap,&newi->cmap);
365: MatSeqDenseSetPreallocation(newi,NULL);
366: if (cpvalues == MAT_COPY_VALUES) {
367: l = (Mat_SeqDense*)newi->data;
368: if (lda>A->rmap->n) {
369: m = A->rmap->n;
370: for (j=0; j<A->cmap->n; j++) {
371: PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
372: }
373: } else {
374: PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
375: }
376: }
377: newi->assembled = PETSC_TRUE;
378: return(0);
379: }
381: static PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
382: {
386: MatCreate(PetscObjectComm((PetscObject)A),newmat);
387: MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
388: MatSetType(*newmat,((PetscObject)A)->type_name);
389: MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
390: return(0);
391: }
394: static PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*);
396: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
397: {
398: MatFactorInfo info;
402: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
403: MatLUFactor_SeqDense(fact,0,0,&info);
404: return(0);
405: }
407: static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
408: {
409: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
410: PetscErrorCode ierr;
411: const PetscScalar *x;
412: PetscScalar *y;
413: PetscBLASInt one = 1,info,m;
416: PetscBLASIntCast(A->rmap->n,&m);
417: VecGetArrayRead(xx,&x);
418: VecGetArray(yy,&y);
419: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
420: if (A->factortype == MAT_FACTOR_LU) {
421: #if defined(PETSC_MISSING_LAPACK_GETRS)
422: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
423: #else
424: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
425: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
426: PetscFPTrapPop();
427: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
428: #endif
429: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
430: #if defined(PETSC_MISSING_LAPACK_POTRS)
431: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
432: #else
433: if (A->spd) {
434: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
435: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
436: PetscFPTrapPop();
437: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
438: #if defined(PETSC_USE_COMPLEX)
439: } else if (A->hermitian) {
440: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
441: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
442: PetscFPTrapPop();
443: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
444: #endif
445: } else { /* symmetric case */
446: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
447: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
448: PetscFPTrapPop();
449: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
450: }
451: #endif
452: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
453: VecRestoreArrayRead(xx,&x);
454: VecRestoreArray(yy,&y);
455: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
456: return(0);
457: }
459: static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
460: {
461: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
463: PetscScalar *b,*x;
464: PetscInt n;
465: PetscBLASInt nrhs,info,m;
466: PetscBool flg;
469: PetscBLASIntCast(A->rmap->n,&m);
470: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
471: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
472: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
473: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
475: MatGetSize(B,NULL,&n);
476: PetscBLASIntCast(n,&nrhs);
477: MatDenseGetArray(B,&b);
478: MatDenseGetArray(X,&x);
480: PetscMemcpy(x,b,m*nrhs*sizeof(PetscScalar));
482: if (A->factortype == MAT_FACTOR_LU) {
483: #if defined(PETSC_MISSING_LAPACK_GETRS)
484: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
485: #else
486: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
487: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
488: PetscFPTrapPop();
489: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
490: #endif
491: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
492: #if defined(PETSC_MISSING_LAPACK_POTRS)
493: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
494: #else
495: if (A->spd) {
496: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
497: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
498: PetscFPTrapPop();
499: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
500: #if defined(PETSC_USE_COMPLEX)
501: } else if (A->hermitian) {
502: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
503: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
504: PetscFPTrapPop();
505: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
506: #endif
507: } else { /* symmetric case */
508: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
509: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
510: PetscFPTrapPop();
511: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
512: }
513: #endif
514: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
516: MatDenseRestoreArray(B,&b);
517: MatDenseRestoreArray(X,&x);
518: PetscLogFlops(nrhs*(2.0*m*m - m));
519: return(0);
520: }
522: static PetscErrorCode MatConjugate_SeqDense(Mat);
524: static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
525: {
526: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
527: PetscErrorCode ierr;
528: const PetscScalar *x;
529: PetscScalar *y;
530: PetscBLASInt one = 1,info,m;
533: PetscBLASIntCast(A->rmap->n,&m);
534: VecGetArrayRead(xx,&x);
535: VecGetArray(yy,&y);
536: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
537: if (A->factortype == MAT_FACTOR_LU) {
538: #if defined(PETSC_MISSING_LAPACK_GETRS)
539: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
540: #else
541: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
542: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
543: PetscFPTrapPop();
544: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
545: #endif
546: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
547: #if defined(PETSC_MISSING_LAPACK_POTRS)
548: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
549: #else
550: if (A->spd) {
551: #if defined(PETSC_USE_COMPLEX)
552: MatConjugate_SeqDense(A);
553: #endif
554: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
555: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
556: PetscFPTrapPop();
557: #if defined(PETSC_USE_COMPLEX)
558: MatConjugate_SeqDense(A);
559: #endif
560: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
561: #if defined(PETSC_USE_COMPLEX)
562: } else if (A->hermitian) {
563: MatConjugate_SeqDense(A);
564: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
565: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
566: PetscFPTrapPop();
567: MatConjugate_SeqDense(A);
568: #endif
569: } else { /* symmetric case */
570: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
571: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
572: PetscFPTrapPop();
573: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
574: }
575: #endif
576: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
577: VecRestoreArrayRead(xx,&x);
578: VecRestoreArray(yy,&y);
579: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
580: return(0);
581: }
583: static PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
584: {
585: PetscErrorCode ierr;
586: const PetscScalar *x;
587: PetscScalar *y,sone = 1.0;
588: Vec tmp = 0;
591: VecGetArrayRead(xx,&x);
592: VecGetArray(yy,&y);
593: if (!A->rmap->n || !A->cmap->n) return(0);
594: if (yy == zz) {
595: VecDuplicate(yy,&tmp);
596: PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
597: VecCopy(yy,tmp);
598: }
599: MatSolve_SeqDense(A,xx,yy);
600: if (tmp) {
601: VecAXPY(yy,sone,tmp);
602: VecDestroy(&tmp);
603: } else {
604: VecAXPY(yy,sone,zz);
605: }
606: VecRestoreArrayRead(xx,&x);
607: VecRestoreArray(yy,&y);
608: PetscLogFlops(A->cmap->n);
609: return(0);
610: }
612: static PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
613: {
614: PetscErrorCode ierr;
615: const PetscScalar *x;
616: PetscScalar *y,sone = 1.0;
617: Vec tmp = NULL;
620: if (!A->rmap->n || !A->cmap->n) return(0);
621: VecGetArrayRead(xx,&x);
622: VecGetArray(yy,&y);
623: if (yy == zz) {
624: VecDuplicate(yy,&tmp);
625: PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
626: VecCopy(yy,tmp);
627: }
628: MatSolveTranspose_SeqDense(A,xx,yy);
629: if (tmp) {
630: VecAXPY(yy,sone,tmp);
631: VecDestroy(&tmp);
632: } else {
633: VecAXPY(yy,sone,zz);
634: }
635: VecRestoreArrayRead(xx,&x);
636: VecRestoreArray(yy,&y);
637: PetscLogFlops(A->cmap->n);
638: return(0);
639: }
641: /* ---------------------------------------------------------------*/
642: /* COMMENT: I have chosen to hide row permutation in the pivots,
643: rather than put it in the Mat->row slot.*/
644: static PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
645: {
646: #if defined(PETSC_MISSING_LAPACK_GETRF)
648: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
649: #else
650: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
652: PetscBLASInt n,m,info;
655: PetscBLASIntCast(A->cmap->n,&n);
656: PetscBLASIntCast(A->rmap->n,&m);
657: if (!mat->pivots) {
658: PetscMalloc1(A->rmap->n,&mat->pivots);
659: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
660: }
661: if (!A->rmap->n || !A->cmap->n) return(0);
662: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
663: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
664: PetscFPTrapPop();
666: if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
667: if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
669: A->ops->solve = MatSolve_SeqDense;
670: A->ops->matsolve = MatMatSolve_SeqDense;
671: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
672: A->ops->solveadd = MatSolveAdd_SeqDense;
673: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
674: A->factortype = MAT_FACTOR_LU;
676: PetscFree(A->solvertype);
677: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
679: PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
680: #endif
681: return(0);
682: }
684: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
685: static PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
686: {
687: #if defined(PETSC_MISSING_LAPACK_POTRF)
689: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
690: #else
691: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
693: PetscBLASInt info,n;
696: PetscBLASIntCast(A->cmap->n,&n);
697: if (!A->rmap->n || !A->cmap->n) return(0);
698: if (A->spd) {
699: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
700: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
701: PetscFPTrapPop();
702: #if defined(PETSC_USE_COMPLEX)
703: } else if (A->hermitian) {
704: if (!mat->pivots) {
705: PetscMalloc1(A->rmap->n,&mat->pivots);
706: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
707: }
708: if (!mat->fwork) {
709: PetscScalar dummy;
711: mat->lfwork = -1;
712: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
713: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
714: PetscFPTrapPop();
715: mat->lfwork = (PetscInt)PetscRealPart(dummy);
716: PetscMalloc1(mat->lfwork,&mat->fwork);
717: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
718: }
719: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
720: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
721: PetscFPTrapPop();
722: #endif
723: } else { /* symmetric case */
724: if (!mat->pivots) {
725: PetscMalloc1(A->rmap->n,&mat->pivots);
726: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
727: }
728: if (!mat->fwork) {
729: PetscScalar dummy;
731: mat->lfwork = -1;
732: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
733: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
734: PetscFPTrapPop();
735: mat->lfwork = (PetscInt)PetscRealPart(dummy);
736: PetscMalloc1(mat->lfwork,&mat->fwork);
737: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
738: }
739: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
740: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
741: PetscFPTrapPop();
742: }
743: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
745: A->ops->solve = MatSolve_SeqDense;
746: A->ops->matsolve = MatMatSolve_SeqDense;
747: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
748: A->ops->solveadd = MatSolveAdd_SeqDense;
749: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
750: A->factortype = MAT_FACTOR_CHOLESKY;
752: PetscFree(A->solvertype);
753: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
755: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
756: #endif
757: return(0);
758: }
761: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
762: {
764: MatFactorInfo info;
767: info.fill = 1.0;
769: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
770: MatCholeskyFactor_SeqDense(fact,0,&info);
771: return(0);
772: }
774: static PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
775: {
777: fact->assembled = PETSC_TRUE;
778: fact->preallocated = PETSC_TRUE;
779: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
780: fact->ops->solve = MatSolve_SeqDense;
781: fact->ops->matsolve = MatMatSolve_SeqDense;
782: fact->ops->solvetranspose = MatSolveTranspose_SeqDense;
783: fact->ops->solveadd = MatSolveAdd_SeqDense;
784: fact->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
785: return(0);
786: }
788: static PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
789: {
791: fact->preallocated = PETSC_TRUE;
792: fact->assembled = PETSC_TRUE;
793: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
794: fact->ops->solve = MatSolve_SeqDense;
795: fact->ops->matsolve = MatMatSolve_SeqDense;
796: fact->ops->solvetranspose = MatSolveTranspose_SeqDense;
797: fact->ops->solveadd = MatSolveAdd_SeqDense;
798: fact->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
799: return(0);
800: }
802: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
803: {
807: MatCreate(PetscObjectComm((PetscObject)A),fact);
808: MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
809: MatSetType(*fact,((PetscObject)A)->type_name);
810: if (ftype == MAT_FACTOR_LU) {
811: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
812: } else {
813: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
814: }
815: (*fact)->factortype = ftype;
817: PetscFree((*fact)->solvertype);
818: PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
819: return(0);
820: }
822: /* ------------------------------------------------------------------*/
823: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
824: {
825: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
826: PetscScalar *x,*v = mat->v,zero = 0.0,xt;
827: const PetscScalar *b;
828: PetscErrorCode ierr;
829: PetscInt m = A->rmap->n,i;
830: PetscBLASInt o = 1,bm;
833: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
834: PetscBLASIntCast(m,&bm);
835: if (flag & SOR_ZERO_INITIAL_GUESS) {
836: /* this is a hack fix, should have another version without the second BLASdotu */
837: VecSet(xx,zero);
838: }
839: VecGetArray(xx,&x);
840: VecGetArrayRead(bb,&b);
841: its = its*lits;
842: if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
843: while (its--) {
844: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
845: for (i=0; i<m; i++) {
846: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
847: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
848: }
849: }
850: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
851: for (i=m-1; i>=0; i--) {
852: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
853: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
854: }
855: }
856: }
857: VecRestoreArrayRead(bb,&b);
858: VecRestoreArray(xx,&x);
859: return(0);
860: }
862: /* -----------------------------------------------------------------*/
863: PETSC_INTERN PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
864: {
865: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
866: const PetscScalar *v = mat->v,*x;
867: PetscScalar *y;
868: PetscErrorCode ierr;
869: PetscBLASInt m, n,_One=1;
870: PetscScalar _DOne=1.0,_DZero=0.0;
873: PetscBLASIntCast(A->rmap->n,&m);
874: PetscBLASIntCast(A->cmap->n,&n);
875: VecGetArrayRead(xx,&x);
876: VecGetArray(yy,&y);
877: if (!A->rmap->n || !A->cmap->n) {
878: PetscBLASInt i;
879: for (i=0; i<n; i++) y[i] = 0.0;
880: } else {
881: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
882: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
883: }
884: VecRestoreArrayRead(xx,&x);
885: VecRestoreArray(yy,&y);
886: return(0);
887: }
889: PETSC_INTERN PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
890: {
891: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
892: PetscScalar *y,_DOne=1.0,_DZero=0.0;
893: PetscErrorCode ierr;
894: PetscBLASInt m, n, _One=1;
895: const PetscScalar *v = mat->v,*x;
898: PetscBLASIntCast(A->rmap->n,&m);
899: PetscBLASIntCast(A->cmap->n,&n);
900: VecGetArrayRead(xx,&x);
901: VecGetArray(yy,&y);
902: if (!A->rmap->n || !A->cmap->n) {
903: PetscBLASInt i;
904: for (i=0; i<m; i++) y[i] = 0.0;
905: } else {
906: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
907: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
908: }
909: VecRestoreArrayRead(xx,&x);
910: VecRestoreArray(yy,&y);
911: return(0);
912: }
914: PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
915: {
916: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
917: const PetscScalar *v = mat->v,*x;
918: PetscScalar *y,_DOne=1.0;
919: PetscErrorCode ierr;
920: PetscBLASInt m, n, _One=1;
923: PetscBLASIntCast(A->rmap->n,&m);
924: PetscBLASIntCast(A->cmap->n,&n);
925: if (!A->rmap->n || !A->cmap->n) return(0);
926: if (zz != yy) {VecCopy(zz,yy);}
927: VecGetArrayRead(xx,&x);
928: VecGetArray(yy,&y);
929: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
930: VecRestoreArrayRead(xx,&x);
931: VecRestoreArray(yy,&y);
932: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
933: return(0);
934: }
936: static PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
937: {
938: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
939: const PetscScalar *v = mat->v,*x;
940: PetscScalar *y;
941: PetscErrorCode ierr;
942: PetscBLASInt m, n, _One=1;
943: PetscScalar _DOne=1.0;
946: PetscBLASIntCast(A->rmap->n,&m);
947: PetscBLASIntCast(A->cmap->n,&n);
948: if (!A->rmap->n || !A->cmap->n) return(0);
949: if (zz != yy) {VecCopy(zz,yy);}
950: VecGetArrayRead(xx,&x);
951: VecGetArray(yy,&y);
952: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
953: VecRestoreArrayRead(xx,&x);
954: VecRestoreArray(yy,&y);
955: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
956: return(0);
957: }
959: /* -----------------------------------------------------------------*/
960: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
961: {
962: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
963: PetscScalar *v;
965: PetscInt i;
968: *ncols = A->cmap->n;
969: if (cols) {
970: PetscMalloc1(A->cmap->n+1,cols);
971: for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
972: }
973: if (vals) {
974: PetscMalloc1(A->cmap->n+1,vals);
975: v = mat->v + row;
976: for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
977: }
978: return(0);
979: }
981: static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
982: {
986: if (cols) {PetscFree(*cols);}
987: if (vals) {PetscFree(*vals); }
988: return(0);
989: }
990: /* ----------------------------------------------------------------*/
991: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
992: {
993: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
994: PetscInt i,j,idx=0;
997: if (!mat->roworiented) {
998: if (addv == INSERT_VALUES) {
999: for (j=0; j<n; j++) {
1000: if (indexn[j] < 0) {idx += m; continue;}
1001: #if defined(PETSC_USE_DEBUG)
1002: 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);
1003: #endif
1004: for (i=0; i<m; i++) {
1005: if (indexm[i] < 0) {idx++; continue;}
1006: #if defined(PETSC_USE_DEBUG)
1007: 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);
1008: #endif
1009: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1010: }
1011: }
1012: } else {
1013: for (j=0; j<n; j++) {
1014: if (indexn[j] < 0) {idx += m; continue;}
1015: #if defined(PETSC_USE_DEBUG)
1016: 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);
1017: #endif
1018: for (i=0; i<m; i++) {
1019: if (indexm[i] < 0) {idx++; continue;}
1020: #if defined(PETSC_USE_DEBUG)
1021: 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);
1022: #endif
1023: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1024: }
1025: }
1026: }
1027: } else {
1028: if (addv == INSERT_VALUES) {
1029: for (i=0; i<m; i++) {
1030: if (indexm[i] < 0) { idx += n; continue;}
1031: #if defined(PETSC_USE_DEBUG)
1032: 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);
1033: #endif
1034: for (j=0; j<n; j++) {
1035: if (indexn[j] < 0) { idx++; continue;}
1036: #if defined(PETSC_USE_DEBUG)
1037: 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);
1038: #endif
1039: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1040: }
1041: }
1042: } else {
1043: for (i=0; i<m; i++) {
1044: if (indexm[i] < 0) { idx += n; continue;}
1045: #if defined(PETSC_USE_DEBUG)
1046: 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);
1047: #endif
1048: for (j=0; j<n; j++) {
1049: if (indexn[j] < 0) { idx++; continue;}
1050: #if defined(PETSC_USE_DEBUG)
1051: 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);
1052: #endif
1053: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1054: }
1055: }
1056: }
1057: }
1058: return(0);
1059: }
1061: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1062: {
1063: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1064: PetscInt i,j;
1067: /* row-oriented output */
1068: for (i=0; i<m; i++) {
1069: if (indexm[i] < 0) {v += n;continue;}
1070: 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);
1071: for (j=0; j<n; j++) {
1072: if (indexn[j] < 0) {v++; continue;}
1073: 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);
1074: *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
1075: }
1076: }
1077: return(0);
1078: }
1080: /* -----------------------------------------------------------------*/
1082: static PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
1083: {
1084: Mat_SeqDense *a;
1086: PetscInt *scols,i,j,nz,header[4];
1087: int fd;
1088: PetscMPIInt size;
1089: PetscInt *rowlengths = 0,M,N,*cols,grows,gcols;
1090: PetscScalar *vals,*svals,*v,*w;
1091: MPI_Comm comm;
1092: PetscBool isbinary;
1095: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1096: if (!isbinary) SETERRQ2(PetscObjectComm((PetscObject)newmat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newmat)->type_name);
1098: /* force binary viewer to load .info file if it has not yet done so */
1099: PetscViewerSetUp(viewer);
1100: PetscObjectGetComm((PetscObject)viewer,&comm);
1101: MPI_Comm_size(comm,&size);
1102: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
1103: PetscViewerBinaryGetDescriptor(viewer,&fd);
1104: PetscBinaryRead(fd,header,4,PETSC_INT);
1105: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
1106: M = header[1]; N = header[2]; nz = header[3];
1108: /* set global size if not set already*/
1109: if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
1110: MatSetSizes(newmat,M,N,M,N);
1111: } else {
1112: /* if sizes and type are already set, check if the vector global sizes are correct */
1113: MatGetSize(newmat,&grows,&gcols);
1114: 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);
1115: }
1116: a = (Mat_SeqDense*)newmat->data;
1117: if (!a->user_alloc) {
1118: MatSeqDenseSetPreallocation(newmat,NULL);
1119: }
1121: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
1122: a = (Mat_SeqDense*)newmat->data;
1123: v = a->v;
1124: /* Allocate some temp space to read in the values and then flip them
1125: from row major to column major */
1126: PetscMalloc1(M*N > 0 ? M*N : 1,&w);
1127: /* read in nonzero values */
1128: PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
1129: /* now flip the values and store them in the matrix*/
1130: for (j=0; j<N; j++) {
1131: for (i=0; i<M; i++) {
1132: *v++ =w[i*N+j];
1133: }
1134: }
1135: PetscFree(w);
1136: } else {
1137: /* read row lengths */
1138: PetscMalloc1(M+1,&rowlengths);
1139: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
1141: a = (Mat_SeqDense*)newmat->data;
1142: v = a->v;
1144: /* read column indices and nonzeros */
1145: PetscMalloc1(nz+1,&scols);
1146: cols = scols;
1147: PetscBinaryRead(fd,cols,nz,PETSC_INT);
1148: PetscMalloc1(nz+1,&svals);
1149: vals = svals;
1150: PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);
1152: /* insert into matrix */
1153: for (i=0; i<M; i++) {
1154: for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
1155: svals += rowlengths[i]; scols += rowlengths[i];
1156: }
1157: PetscFree(vals);
1158: PetscFree(cols);
1159: PetscFree(rowlengths);
1160: }
1161: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1162: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1163: return(0);
1164: }
1166: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1167: {
1168: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1169: PetscErrorCode ierr;
1170: PetscInt i,j;
1171: const char *name;
1172: PetscScalar *v;
1173: PetscViewerFormat format;
1174: #if defined(PETSC_USE_COMPLEX)
1175: PetscBool allreal = PETSC_TRUE;
1176: #endif
1179: PetscViewerGetFormat(viewer,&format);
1180: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1181: return(0); /* do nothing for now */
1182: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1183: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1184: for (i=0; i<A->rmap->n; i++) {
1185: v = a->v + i;
1186: PetscViewerASCIIPrintf(viewer,"row %D:",i);
1187: for (j=0; j<A->cmap->n; j++) {
1188: #if defined(PETSC_USE_COMPLEX)
1189: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1190: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1191: } else if (PetscRealPart(*v)) {
1192: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1193: }
1194: #else
1195: if (*v) {
1196: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1197: }
1198: #endif
1199: v += a->lda;
1200: }
1201: PetscViewerASCIIPrintf(viewer,"\n");
1202: }
1203: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1204: } else {
1205: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1206: #if defined(PETSC_USE_COMPLEX)
1207: /* determine if matrix has all real values */
1208: v = a->v;
1209: for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1210: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1211: }
1212: #endif
1213: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1214: PetscObjectGetName((PetscObject)A,&name);
1215: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1216: PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1217: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1218: }
1220: for (i=0; i<A->rmap->n; i++) {
1221: v = a->v + i;
1222: for (j=0; j<A->cmap->n; j++) {
1223: #if defined(PETSC_USE_COMPLEX)
1224: if (allreal) {
1225: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1226: } else {
1227: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1228: }
1229: #else
1230: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1231: #endif
1232: v += a->lda;
1233: }
1234: PetscViewerASCIIPrintf(viewer,"\n");
1235: }
1236: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1237: PetscViewerASCIIPrintf(viewer,"];\n");
1238: }
1239: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1240: }
1241: PetscViewerFlush(viewer);
1242: return(0);
1243: }
1245: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1246: {
1247: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1248: PetscErrorCode ierr;
1249: int fd;
1250: PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1251: PetscScalar *v,*anonz,*vals;
1252: PetscViewerFormat format;
1255: PetscViewerBinaryGetDescriptor(viewer,&fd);
1257: PetscViewerGetFormat(viewer,&format);
1258: if (format == PETSC_VIEWER_NATIVE) {
1259: /* store the matrix as a dense matrix */
1260: PetscMalloc1(4,&col_lens);
1262: col_lens[0] = MAT_FILE_CLASSID;
1263: col_lens[1] = m;
1264: col_lens[2] = n;
1265: col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;
1267: PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1268: PetscFree(col_lens);
1270: /* write out matrix, by rows */
1271: PetscMalloc1(m*n+1,&vals);
1272: v = a->v;
1273: for (j=0; j<n; j++) {
1274: for (i=0; i<m; i++) {
1275: vals[j + i*n] = *v++;
1276: }
1277: }
1278: PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1279: PetscFree(vals);
1280: } else {
1281: PetscMalloc1(4+nz,&col_lens);
1283: col_lens[0] = MAT_FILE_CLASSID;
1284: col_lens[1] = m;
1285: col_lens[2] = n;
1286: col_lens[3] = nz;
1288: /* store lengths of each row and write (including header) to file */
1289: for (i=0; i<m; i++) col_lens[4+i] = n;
1290: PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);
1292: /* Possibly should write in smaller increments, not whole matrix at once? */
1293: /* store column indices (zero start index) */
1294: ict = 0;
1295: for (i=0; i<m; i++) {
1296: for (j=0; j<n; j++) col_lens[ict++] = j;
1297: }
1298: PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1299: PetscFree(col_lens);
1301: /* store nonzero values */
1302: PetscMalloc1(nz+1,&anonz);
1303: ict = 0;
1304: for (i=0; i<m; i++) {
1305: v = a->v + i;
1306: for (j=0; j<n; j++) {
1307: anonz[ict++] = *v; v += a->lda;
1308: }
1309: }
1310: PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1311: PetscFree(anonz);
1312: }
1313: return(0);
1314: }
1316: #include <petscdraw.h>
1317: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1318: {
1319: Mat A = (Mat) Aa;
1320: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1321: PetscErrorCode ierr;
1322: PetscInt m = A->rmap->n,n = A->cmap->n,i,j;
1323: int color = PETSC_DRAW_WHITE;
1324: PetscScalar *v = a->v;
1325: PetscViewer viewer;
1326: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1327: PetscViewerFormat format;
1330: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1331: PetscViewerGetFormat(viewer,&format);
1332: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
1334: /* Loop over matrix elements drawing boxes */
1336: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1337: PetscDrawCollectiveBegin(draw);
1338: /* Blue for negative and Red for positive */
1339: for (j = 0; j < n; j++) {
1340: x_l = j; x_r = x_l + 1.0;
1341: for (i = 0; i < m; i++) {
1342: y_l = m - i - 1.0;
1343: y_r = y_l + 1.0;
1344: if (PetscRealPart(v[j*m+i]) > 0.) {
1345: color = PETSC_DRAW_RED;
1346: } else if (PetscRealPart(v[j*m+i]) < 0.) {
1347: color = PETSC_DRAW_BLUE;
1348: } else {
1349: continue;
1350: }
1351: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1352: }
1353: }
1354: PetscDrawCollectiveEnd(draw);
1355: } else {
1356: /* use contour shading to indicate magnitude of values */
1357: /* first determine max of all nonzero values */
1358: PetscReal minv = 0.0, maxv = 0.0;
1359: PetscDraw popup;
1361: for (i=0; i < m*n; i++) {
1362: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1363: }
1364: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1365: PetscDrawGetPopup(draw,&popup);
1366: PetscDrawScalePopup(popup,minv,maxv);
1368: PetscDrawCollectiveBegin(draw);
1369: for (j=0; j<n; j++) {
1370: x_l = j;
1371: x_r = x_l + 1.0;
1372: for (i=0; i<m; i++) {
1373: y_l = m - i - 1.0;
1374: y_r = y_l + 1.0;
1375: color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1376: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1377: }
1378: }
1379: PetscDrawCollectiveEnd(draw);
1380: }
1381: return(0);
1382: }
1384: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1385: {
1386: PetscDraw draw;
1387: PetscBool isnull;
1388: PetscReal xr,yr,xl,yl,h,w;
1392: PetscViewerDrawGetDraw(viewer,0,&draw);
1393: PetscDrawIsNull(draw,&isnull);
1394: if (isnull) return(0);
1396: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1397: xr += w; yr += h; xl = -w; yl = -h;
1398: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1399: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1400: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1401: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1402: PetscDrawSave(draw);
1403: return(0);
1404: }
1406: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1407: {
1409: PetscBool iascii,isbinary,isdraw;
1412: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1413: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1414: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1416: if (iascii) {
1417: MatView_SeqDense_ASCII(A,viewer);
1418: } else if (isbinary) {
1419: MatView_SeqDense_Binary(A,viewer);
1420: } else if (isdraw) {
1421: MatView_SeqDense_Draw(A,viewer);
1422: }
1423: return(0);
1424: }
1426: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar array[])
1427: {
1428: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1431: a->unplacedarray = a->v;
1432: a->unplaced_user_alloc = a->user_alloc;
1433: a->v = (PetscScalar*) array;
1434: return(0);
1435: }
1437: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1438: {
1439: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1442: a->v = a->unplacedarray;
1443: a->user_alloc = a->unplaced_user_alloc;
1444: a->unplacedarray = NULL;
1445: return(0);
1446: }
1448: static PetscErrorCode MatDestroy_SeqDense(Mat mat)
1449: {
1450: Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
1454: #if defined(PETSC_USE_LOG)
1455: PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1456: #endif
1457: PetscFree(l->pivots);
1458: PetscFree(l->fwork);
1459: MatDestroy(&l->ptapwork);
1460: if (!l->user_alloc) {PetscFree(l->v);}
1461: PetscFree(mat->data);
1463: PetscObjectChangeTypeName((PetscObject)mat,0);
1464: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1465: PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1466: PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1467: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1468: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1469: #if defined(PETSC_HAVE_ELEMENTAL)
1470: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1471: #endif
1472: PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1473: PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1474: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1475: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1476: PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);
1477: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1478: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1479: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1480: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
1481: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
1482: return(0);
1483: }
1485: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1486: {
1487: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1489: PetscInt k,j,m,n,M;
1490: PetscScalar *v,tmp;
1493: v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1494: if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */
1495: for (j=0; j<m; j++) {
1496: for (k=0; k<j; k++) {
1497: tmp = v[j + k*M];
1498: v[j + k*M] = v[k + j*M];
1499: v[k + j*M] = tmp;
1500: }
1501: }
1502: } else { /* out-of-place transpose */
1503: Mat tmat;
1504: Mat_SeqDense *tmatd;
1505: PetscScalar *v2;
1506: PetscInt M2;
1508: if (reuse != MAT_REUSE_MATRIX) {
1509: MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1510: MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1511: MatSetType(tmat,((PetscObject)A)->type_name);
1512: MatSeqDenseSetPreallocation(tmat,NULL);
1513: } else {
1514: tmat = *matout;
1515: }
1516: tmatd = (Mat_SeqDense*)tmat->data;
1517: v = mat->v; v2 = tmatd->v; M2 = tmatd->lda;
1518: for (j=0; j<n; j++) {
1519: for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1520: }
1521: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1522: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1523: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat;
1524: else {
1525: MatHeaderMerge(A,&tmat);
1526: }
1527: }
1528: return(0);
1529: }
1531: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg)
1532: {
1533: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1534: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1535: PetscInt i,j;
1536: PetscScalar *v1,*v2;
1539: if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1540: if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1541: for (i=0; i<A1->rmap->n; i++) {
1542: v1 = mat1->v+i; v2 = mat2->v+i;
1543: for (j=0; j<A1->cmap->n; j++) {
1544: if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1545: v1 += mat1->lda; v2 += mat2->lda;
1546: }
1547: }
1548: *flg = PETSC_TRUE;
1549: return(0);
1550: }
1552: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1553: {
1554: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1556: PetscInt i,n,len;
1557: PetscScalar *x,zero = 0.0;
1560: VecSet(v,zero);
1561: VecGetSize(v,&n);
1562: VecGetArray(v,&x);
1563: len = PetscMin(A->rmap->n,A->cmap->n);
1564: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1565: for (i=0; i<len; i++) {
1566: x[i] = mat->v[i*mat->lda + i];
1567: }
1568: VecRestoreArray(v,&x);
1569: return(0);
1570: }
1572: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1573: {
1574: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1575: const PetscScalar *l,*r;
1576: PetscScalar x,*v;
1577: PetscErrorCode ierr;
1578: PetscInt i,j,m = A->rmap->n,n = A->cmap->n;
1581: if (ll) {
1582: VecGetSize(ll,&m);
1583: VecGetArrayRead(ll,&l);
1584: if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1585: for (i=0; i<m; i++) {
1586: x = l[i];
1587: v = mat->v + i;
1588: for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1589: }
1590: VecRestoreArrayRead(ll,&l);
1591: PetscLogFlops(1.0*n*m);
1592: }
1593: if (rr) {
1594: VecGetSize(rr,&n);
1595: VecGetArrayRead(rr,&r);
1596: if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1597: for (i=0; i<n; i++) {
1598: x = r[i];
1599: v = mat->v + i*mat->lda;
1600: for (j=0; j<m; j++) (*v++) *= x;
1601: }
1602: VecRestoreArrayRead(rr,&r);
1603: PetscLogFlops(1.0*n*m);
1604: }
1605: return(0);
1606: }
1608: static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1609: {
1610: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1611: PetscScalar *v = mat->v;
1612: PetscReal sum = 0.0;
1613: PetscInt lda =mat->lda,m=A->rmap->n,i,j;
1617: if (type == NORM_FROBENIUS) {
1618: if (lda>m) {
1619: for (j=0; j<A->cmap->n; j++) {
1620: v = mat->v+j*lda;
1621: for (i=0; i<m; i++) {
1622: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1623: }
1624: }
1625: } else {
1626: #if defined(PETSC_USE_REAL___FP16)
1627: PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1628: *nrm = BLASnrm2_(&cnt,v,&one);
1629: }
1630: #else
1631: for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1632: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1633: }
1634: }
1635: *nrm = PetscSqrtReal(sum);
1636: #endif
1637: PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1638: } else if (type == NORM_1) {
1639: *nrm = 0.0;
1640: for (j=0; j<A->cmap->n; j++) {
1641: v = mat->v + j*mat->lda;
1642: sum = 0.0;
1643: for (i=0; i<A->rmap->n; i++) {
1644: sum += PetscAbsScalar(*v); v++;
1645: }
1646: if (sum > *nrm) *nrm = sum;
1647: }
1648: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1649: } else if (type == NORM_INFINITY) {
1650: *nrm = 0.0;
1651: for (j=0; j<A->rmap->n; j++) {
1652: v = mat->v + j;
1653: sum = 0.0;
1654: for (i=0; i<A->cmap->n; i++) {
1655: sum += PetscAbsScalar(*v); v += mat->lda;
1656: }
1657: if (sum > *nrm) *nrm = sum;
1658: }
1659: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1660: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1661: return(0);
1662: }
1664: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1665: {
1666: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
1670: switch (op) {
1671: case MAT_ROW_ORIENTED:
1672: aij->roworiented = flg;
1673: break;
1674: case MAT_NEW_NONZERO_LOCATIONS:
1675: case MAT_NEW_NONZERO_LOCATION_ERR:
1676: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1677: case MAT_NEW_DIAGONALS:
1678: case MAT_KEEP_NONZERO_PATTERN:
1679: case MAT_IGNORE_OFF_PROC_ENTRIES:
1680: case MAT_USE_HASH_TABLE:
1681: case MAT_IGNORE_ZERO_ENTRIES:
1682: case MAT_IGNORE_LOWER_TRIANGULAR:
1683: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1684: break;
1685: case MAT_SPD:
1686: case MAT_SYMMETRIC:
1687: case MAT_STRUCTURALLY_SYMMETRIC:
1688: case MAT_HERMITIAN:
1689: case MAT_SYMMETRY_ETERNAL:
1690: /* These options are handled directly by MatSetOption() */
1691: break;
1692: default:
1693: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1694: }
1695: return(0);
1696: }
1698: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1699: {
1700: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1702: PetscInt lda=l->lda,m=A->rmap->n,j;
1705: if (lda>m) {
1706: for (j=0; j<A->cmap->n; j++) {
1707: PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1708: }
1709: } else {
1710: PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1711: }
1712: return(0);
1713: }
1715: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1716: {
1717: PetscErrorCode ierr;
1718: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1719: PetscInt m = l->lda, n = A->cmap->n, i,j;
1720: PetscScalar *slot,*bb;
1721: const PetscScalar *xx;
1724: #if defined(PETSC_USE_DEBUG)
1725: for (i=0; i<N; i++) {
1726: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1727: 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);
1728: }
1729: #endif
1731: /* fix right hand side if needed */
1732: if (x && b) {
1733: VecGetArrayRead(x,&xx);
1734: VecGetArray(b,&bb);
1735: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1736: VecRestoreArrayRead(x,&xx);
1737: VecRestoreArray(b,&bb);
1738: }
1740: for (i=0; i<N; i++) {
1741: slot = l->v + rows[i];
1742: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1743: }
1744: if (diag != 0.0) {
1745: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1746: for (i=0; i<N; i++) {
1747: slot = l->v + (m+1)*rows[i];
1748: *slot = diag;
1749: }
1750: }
1751: return(0);
1752: }
1754: static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1755: {
1756: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1759: 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");
1760: *array = mat->v;
1761: return(0);
1762: }
1764: static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1765: {
1767: return(0);
1768: }
1770: /*@C
1771: MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored
1773: Logically Collective on Mat
1775: Input Parameter:
1776: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1778: Output Parameter:
1779: . array - pointer to the data
1781: Level: intermediate
1783: .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
1784: @*/
1785: PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array)
1786: {
1790: PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1791: return(0);
1792: }
1794: /*@C
1795: MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
1797: Logically Collective on Mat
1799: Input Parameters:
1800: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1801: . array - pointer to the data
1803: Level: intermediate
1805: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
1806: @*/
1807: PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array)
1808: {
1812: PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1813: if (array) *array = NULL;
1814: PetscObjectStateIncrease((PetscObject)A);
1815: return(0);
1816: }
1818: /*@C
1819: MatDenseGetArrayRead - gives access to the array where the data for a SeqDense matrix is stored
1821: Not Collective
1823: Input Parameter:
1824: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1826: Output Parameter:
1827: . array - pointer to the data
1829: Level: intermediate
1831: .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead()
1832: @*/
1833: PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array)
1834: {
1838: PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
1839: return(0);
1840: }
1842: /*@C
1843: MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
1845: Not Collective
1847: Input Parameters:
1848: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1849: . array - pointer to the data
1851: Level: intermediate
1853: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray()
1854: @*/
1855: PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
1856: {
1860: PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
1861: if (array) *array = NULL;
1862: return(0);
1863: }
1865: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1866: {
1867: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1869: PetscInt i,j,nrows,ncols;
1870: const PetscInt *irow,*icol;
1871: PetscScalar *av,*bv,*v = mat->v;
1872: Mat newmat;
1875: ISGetIndices(isrow,&irow);
1876: ISGetIndices(iscol,&icol);
1877: ISGetLocalSize(isrow,&nrows);
1878: ISGetLocalSize(iscol,&ncols);
1880: /* Check submatrixcall */
1881: if (scall == MAT_REUSE_MATRIX) {
1882: PetscInt n_cols,n_rows;
1883: MatGetSize(*B,&n_rows,&n_cols);
1884: if (n_rows != nrows || n_cols != ncols) {
1885: /* resize the result matrix to match number of requested rows/columns */
1886: MatSetSizes(*B,nrows,ncols,nrows,ncols);
1887: }
1888: newmat = *B;
1889: } else {
1890: /* Create and fill new matrix */
1891: MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1892: MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1893: MatSetType(newmat,((PetscObject)A)->type_name);
1894: MatSeqDenseSetPreallocation(newmat,NULL);
1895: }
1897: /* Now extract the data pointers and do the copy,column at a time */
1898: bv = ((Mat_SeqDense*)newmat->data)->v;
1900: for (i=0; i<ncols; i++) {
1901: av = v + mat->lda*icol[i];
1902: for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1903: }
1905: /* Assemble the matrices so that the correct flags are set */
1906: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1907: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1909: /* Free work space */
1910: ISRestoreIndices(isrow,&irow);
1911: ISRestoreIndices(iscol,&icol);
1912: *B = newmat;
1913: return(0);
1914: }
1916: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
1917: {
1919: PetscInt i;
1922: if (scall == MAT_INITIAL_MATRIX) {
1923: PetscCalloc1(n+1,B);
1924: }
1926: for (i=0; i<n; i++) {
1927: MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1928: }
1929: return(0);
1930: }
1932: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1933: {
1935: return(0);
1936: }
1938: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1939: {
1941: return(0);
1942: }
1944: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1945: {
1946: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1948: PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;
1951: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1952: if (A->ops->copy != B->ops->copy) {
1953: MatCopy_Basic(A,B,str);
1954: return(0);
1955: }
1956: if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1957: if (lda1>m || lda2>m) {
1958: for (j=0; j<n; j++) {
1959: PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1960: }
1961: } else {
1962: PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1963: }
1964: PetscObjectStateIncrease((PetscObject)B);
1965: return(0);
1966: }
1968: static PetscErrorCode MatSetUp_SeqDense(Mat A)
1969: {
1973: MatSeqDenseSetPreallocation(A,0);
1974: return(0);
1975: }
1977: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1978: {
1979: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1980: PetscInt i,nz = A->rmap->n*A->cmap->n;
1981: PetscScalar *aa = a->v;
1984: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1985: return(0);
1986: }
1988: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1989: {
1990: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1991: PetscInt i,nz = A->rmap->n*A->cmap->n;
1992: PetscScalar *aa = a->v;
1995: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1996: return(0);
1997: }
1999: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2000: {
2001: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2002: PetscInt i,nz = A->rmap->n*A->cmap->n;
2003: PetscScalar *aa = a->v;
2006: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2007: return(0);
2008: }
2010: /* ----------------------------------------------------------------*/
2011: PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2012: {
2016: if (scall == MAT_INITIAL_MATRIX) {
2017: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
2018: MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2019: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
2020: }
2021: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
2022: MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2023: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
2024: return(0);
2025: }
2027: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2028: {
2030: PetscInt m=A->rmap->n,n=B->cmap->n;
2031: Mat Cmat;
2034: 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);
2035: MatCreate(PETSC_COMM_SELF,&Cmat);
2036: MatSetSizes(Cmat,m,n,m,n);
2037: MatSetType(Cmat,MATSEQDENSE);
2038: MatSeqDenseSetPreallocation(Cmat,NULL);
2040: *C = Cmat;
2041: return(0);
2042: }
2044: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2045: {
2046: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2047: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2048: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2049: PetscBLASInt m,n,k;
2050: PetscScalar _DOne=1.0,_DZero=0.0;
2052: PetscBool flg;
2055: PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
2056: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Second matrix must be dense");
2058: /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/
2059: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
2060: if (flg) {
2061: C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
2062: (*C->ops->matmultnumeric)(A,B,C);
2063: return(0);
2064: }
2066: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
2067: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense");
2068: PetscBLASIntCast(C->rmap->n,&m);
2069: PetscBLASIntCast(C->cmap->n,&n);
2070: PetscBLASIntCast(A->cmap->n,&k);
2071: if (!m || !n || !k) return(0);
2072: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2073: return(0);
2074: }
2076: PetscErrorCode MatMatTransposeMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2077: {
2081: if (scall == MAT_INITIAL_MATRIX) {
2082: PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);
2083: MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2084: PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);
2085: }
2086: PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);
2087: MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);
2088: PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);
2089: return(0);
2090: }
2092: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2093: {
2095: PetscInt m=A->rmap->n,n=B->rmap->n;
2096: Mat Cmat;
2099: 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);
2100: MatCreate(PETSC_COMM_SELF,&Cmat);
2101: MatSetSizes(Cmat,m,n,m,n);
2102: MatSetType(Cmat,MATSEQDENSE);
2103: MatSeqDenseSetPreallocation(Cmat,NULL);
2105: Cmat->assembled = PETSC_TRUE;
2107: *C = Cmat;
2108: return(0);
2109: }
2111: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2112: {
2113: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2114: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2115: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2116: PetscBLASInt m,n,k;
2117: PetscScalar _DOne=1.0,_DZero=0.0;
2121: PetscBLASIntCast(C->rmap->n,&m);
2122: PetscBLASIntCast(C->cmap->n,&n);
2123: PetscBLASIntCast(A->cmap->n,&k);
2124: if (!m || !n || !k) return(0);
2125: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2126: return(0);
2127: }
2129: PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2130: {
2134: if (scall == MAT_INITIAL_MATRIX) {
2135: PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
2136: MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2137: PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
2138: }
2139: PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
2140: MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2141: PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
2142: return(0);
2143: }
2145: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2146: {
2148: PetscInt m=A->cmap->n,n=B->cmap->n;
2149: Mat Cmat;
2152: 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);
2153: MatCreate(PETSC_COMM_SELF,&Cmat);
2154: MatSetSizes(Cmat,m,n,m,n);
2155: MatSetType(Cmat,MATSEQDENSE);
2156: MatSeqDenseSetPreallocation(Cmat,NULL);
2158: Cmat->assembled = PETSC_TRUE;
2160: *C = Cmat;
2161: return(0);
2162: }
2164: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2165: {
2166: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2167: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2168: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2169: PetscBLASInt m,n,k;
2170: PetscScalar _DOne=1.0,_DZero=0.0;
2174: PetscBLASIntCast(C->rmap->n,&m);
2175: PetscBLASIntCast(C->cmap->n,&n);
2176: PetscBLASIntCast(A->rmap->n,&k);
2177: if (!m || !n || !k) return(0);
2178: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2179: return(0);
2180: }
2182: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2183: {
2184: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2186: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2187: PetscScalar *x;
2188: MatScalar *aa = a->v;
2191: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2193: VecSet(v,0.0);
2194: VecGetArray(v,&x);
2195: VecGetLocalSize(v,&p);
2196: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2197: for (i=0; i<m; i++) {
2198: x[i] = aa[i]; if (idx) idx[i] = 0;
2199: for (j=1; j<n; j++) {
2200: if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2201: }
2202: }
2203: VecRestoreArray(v,&x);
2204: return(0);
2205: }
2207: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2208: {
2209: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2211: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2212: PetscScalar *x;
2213: PetscReal atmp;
2214: MatScalar *aa = a->v;
2217: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2219: VecSet(v,0.0);
2220: VecGetArray(v,&x);
2221: VecGetLocalSize(v,&p);
2222: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2223: for (i=0; i<m; i++) {
2224: x[i] = PetscAbsScalar(aa[i]);
2225: for (j=1; j<n; j++) {
2226: atmp = PetscAbsScalar(aa[i+m*j]);
2227: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2228: }
2229: }
2230: VecRestoreArray(v,&x);
2231: return(0);
2232: }
2234: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2235: {
2236: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2238: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2239: PetscScalar *x;
2240: MatScalar *aa = a->v;
2243: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2245: VecSet(v,0.0);
2246: VecGetArray(v,&x);
2247: VecGetLocalSize(v,&p);
2248: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2249: for (i=0; i<m; i++) {
2250: x[i] = aa[i]; if (idx) idx[i] = 0;
2251: for (j=1; j<n; j++) {
2252: if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2253: }
2254: }
2255: VecRestoreArray(v,&x);
2256: return(0);
2257: }
2259: static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2260: {
2261: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2263: PetscScalar *x;
2266: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2268: VecGetArray(v,&x);
2269: PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
2270: VecRestoreArray(v,&x);
2271: return(0);
2272: }
2274: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2275: {
2277: PetscInt i,j,m,n;
2278: PetscScalar *a;
2281: MatGetSize(A,&m,&n);
2282: PetscMemzero(norms,n*sizeof(PetscReal));
2283: MatDenseGetArray(A,&a);
2284: if (type == NORM_2) {
2285: for (i=0; i<n; i++) {
2286: for (j=0; j<m; j++) {
2287: norms[i] += PetscAbsScalar(a[j]*a[j]);
2288: }
2289: a += m;
2290: }
2291: } else if (type == NORM_1) {
2292: for (i=0; i<n; i++) {
2293: for (j=0; j<m; j++) {
2294: norms[i] += PetscAbsScalar(a[j]);
2295: }
2296: a += m;
2297: }
2298: } else if (type == NORM_INFINITY) {
2299: for (i=0; i<n; i++) {
2300: for (j=0; j<m; j++) {
2301: norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2302: }
2303: a += m;
2304: }
2305: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2306: MatDenseRestoreArray(A,&a);
2307: if (type == NORM_2) {
2308: for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2309: }
2310: return(0);
2311: }
2313: static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2314: {
2316: PetscScalar *a;
2317: PetscInt m,n,i;
2320: MatGetSize(x,&m,&n);
2321: MatDenseGetArray(x,&a);
2322: for (i=0; i<m*n; i++) {
2323: PetscRandomGetValue(rctx,a+i);
2324: }
2325: MatDenseRestoreArray(x,&a);
2326: return(0);
2327: }
2329: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d)
2330: {
2332: *missing = PETSC_FALSE;
2333: return(0);
2334: }
2336: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2337: {
2338: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2341: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2342: *vals = a->v+col*a->lda;
2343: return(0);
2344: }
2346: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2347: {
2349: *vals = 0; /* user cannot accidently use the array later */
2350: return(0);
2351: }
2353: /* -------------------------------------------------------------------*/
2354: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2355: MatGetRow_SeqDense,
2356: MatRestoreRow_SeqDense,
2357: MatMult_SeqDense,
2358: /* 4*/ MatMultAdd_SeqDense,
2359: MatMultTranspose_SeqDense,
2360: MatMultTransposeAdd_SeqDense,
2361: 0,
2362: 0,
2363: 0,
2364: /* 10*/ 0,
2365: MatLUFactor_SeqDense,
2366: MatCholeskyFactor_SeqDense,
2367: MatSOR_SeqDense,
2368: MatTranspose_SeqDense,
2369: /* 15*/ MatGetInfo_SeqDense,
2370: MatEqual_SeqDense,
2371: MatGetDiagonal_SeqDense,
2372: MatDiagonalScale_SeqDense,
2373: MatNorm_SeqDense,
2374: /* 20*/ MatAssemblyBegin_SeqDense,
2375: MatAssemblyEnd_SeqDense,
2376: MatSetOption_SeqDense,
2377: MatZeroEntries_SeqDense,
2378: /* 24*/ MatZeroRows_SeqDense,
2379: 0,
2380: 0,
2381: 0,
2382: 0,
2383: /* 29*/ MatSetUp_SeqDense,
2384: 0,
2385: 0,
2386: 0,
2387: 0,
2388: /* 34*/ MatDuplicate_SeqDense,
2389: 0,
2390: 0,
2391: 0,
2392: 0,
2393: /* 39*/ MatAXPY_SeqDense,
2394: MatCreateSubMatrices_SeqDense,
2395: 0,
2396: MatGetValues_SeqDense,
2397: MatCopy_SeqDense,
2398: /* 44*/ MatGetRowMax_SeqDense,
2399: MatScale_SeqDense,
2400: MatShift_Basic,
2401: 0,
2402: MatZeroRowsColumns_SeqDense,
2403: /* 49*/ MatSetRandom_SeqDense,
2404: 0,
2405: 0,
2406: 0,
2407: 0,
2408: /* 54*/ 0,
2409: 0,
2410: 0,
2411: 0,
2412: 0,
2413: /* 59*/ 0,
2414: MatDestroy_SeqDense,
2415: MatView_SeqDense,
2416: 0,
2417: 0,
2418: /* 64*/ 0,
2419: 0,
2420: 0,
2421: 0,
2422: 0,
2423: /* 69*/ MatGetRowMaxAbs_SeqDense,
2424: 0,
2425: 0,
2426: 0,
2427: 0,
2428: /* 74*/ 0,
2429: 0,
2430: 0,
2431: 0,
2432: 0,
2433: /* 79*/ 0,
2434: 0,
2435: 0,
2436: 0,
2437: /* 83*/ MatLoad_SeqDense,
2438: 0,
2439: MatIsHermitian_SeqDense,
2440: 0,
2441: 0,
2442: 0,
2443: /* 89*/ MatMatMult_SeqDense_SeqDense,
2444: MatMatMultSymbolic_SeqDense_SeqDense,
2445: MatMatMultNumeric_SeqDense_SeqDense,
2446: MatPtAP_SeqDense_SeqDense,
2447: MatPtAPSymbolic_SeqDense_SeqDense,
2448: /* 94*/ MatPtAPNumeric_SeqDense_SeqDense,
2449: MatMatTransposeMult_SeqDense_SeqDense,
2450: MatMatTransposeMultSymbolic_SeqDense_SeqDense,
2451: MatMatTransposeMultNumeric_SeqDense_SeqDense,
2452: 0,
2453: /* 99*/ 0,
2454: 0,
2455: 0,
2456: MatConjugate_SeqDense,
2457: 0,
2458: /*104*/ 0,
2459: MatRealPart_SeqDense,
2460: MatImaginaryPart_SeqDense,
2461: 0,
2462: 0,
2463: /*109*/ 0,
2464: 0,
2465: MatGetRowMin_SeqDense,
2466: MatGetColumnVector_SeqDense,
2467: MatMissingDiagonal_SeqDense,
2468: /*114*/ 0,
2469: 0,
2470: 0,
2471: 0,
2472: 0,
2473: /*119*/ 0,
2474: 0,
2475: 0,
2476: 0,
2477: 0,
2478: /*124*/ 0,
2479: MatGetColumnNorms_SeqDense,
2480: 0,
2481: 0,
2482: 0,
2483: /*129*/ 0,
2484: MatTransposeMatMult_SeqDense_SeqDense,
2485: MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2486: MatTransposeMatMultNumeric_SeqDense_SeqDense,
2487: 0,
2488: /*134*/ 0,
2489: 0,
2490: 0,
2491: 0,
2492: 0,
2493: /*139*/ 0,
2494: 0,
2495: 0,
2496: 0,
2497: 0,
2498: /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqDense
2499: };
2501: /*@C
2502: MatCreateSeqDense - Creates a sequential dense matrix that
2503: is stored in column major order (the usual Fortran 77 manner). Many
2504: of the matrix operations use the BLAS and LAPACK routines.
2506: Collective on MPI_Comm
2508: Input Parameters:
2509: + comm - MPI communicator, set to PETSC_COMM_SELF
2510: . m - number of rows
2511: . n - number of columns
2512: - data - optional location of matrix data in column major order. Set data=NULL for PETSc
2513: to control all matrix memory allocation.
2515: Output Parameter:
2516: . A - the matrix
2518: Notes:
2519: The data input variable is intended primarily for Fortran programmers
2520: who wish to allocate their own matrix memory space. Most users should
2521: set data=NULL.
2523: Level: intermediate
2525: .keywords: dense, matrix, LAPACK, BLAS
2527: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2528: @*/
2529: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2530: {
2534: MatCreate(comm,A);
2535: MatSetSizes(*A,m,n,m,n);
2536: MatSetType(*A,MATSEQDENSE);
2537: MatSeqDenseSetPreallocation(*A,data);
2538: return(0);
2539: }
2541: /*@C
2542: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
2544: Collective on MPI_Comm
2546: Input Parameters:
2547: + B - the matrix
2548: - data - the array (or NULL)
2550: Notes:
2551: The data input variable is intended primarily for Fortran programmers
2552: who wish to allocate their own matrix memory space. Most users should
2553: need not call this routine.
2555: Level: intermediate
2557: .keywords: dense, matrix, LAPACK, BLAS
2559: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA()
2561: @*/
2562: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2563: {
2567: PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2568: return(0);
2569: }
2571: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2572: {
2573: Mat_SeqDense *b;
2577: B->preallocated = PETSC_TRUE;
2579: PetscLayoutSetUp(B->rmap);
2580: PetscLayoutSetUp(B->cmap);
2582: b = (Mat_SeqDense*)B->data;
2583: b->Mmax = B->rmap->n;
2584: b->Nmax = B->cmap->n;
2585: if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;
2587: PetscIntMultError(b->lda,b->Nmax,NULL);
2588: if (!data) { /* petsc-allocated storage */
2589: if (!b->user_alloc) { PetscFree(b->v); }
2590: PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2591: PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));
2593: b->user_alloc = PETSC_FALSE;
2594: } else { /* user-allocated storage */
2595: if (!b->user_alloc) { PetscFree(b->v); }
2596: b->v = data;
2597: b->user_alloc = PETSC_TRUE;
2598: }
2599: B->assembled = PETSC_TRUE;
2600: return(0);
2601: }
2603: #if defined(PETSC_HAVE_ELEMENTAL)
2604: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2605: {
2606: Mat mat_elemental;
2608: PetscScalar *array,*v_colwise;
2609: PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;
2612: PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2613: MatDenseGetArray(A,&array);
2614: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2615: k = 0;
2616: for (j=0; j<N; j++) {
2617: cols[j] = j;
2618: for (i=0; i<M; i++) {
2619: v_colwise[j*M+i] = array[k++];
2620: }
2621: }
2622: for (i=0; i<M; i++) {
2623: rows[i] = i;
2624: }
2625: MatDenseRestoreArray(A,&array);
2627: MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2628: MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2629: MatSetType(mat_elemental,MATELEMENTAL);
2630: MatSetUp(mat_elemental);
2632: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2633: MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2634: MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2635: MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2636: PetscFree3(v_colwise,rows,cols);
2638: if (reuse == MAT_INPLACE_MATRIX) {
2639: MatHeaderReplace(A,&mat_elemental);
2640: } else {
2641: *newmat = mat_elemental;
2642: }
2643: return(0);
2644: }
2645: #endif
2647: /*@C
2648: MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array
2650: Input parameter:
2651: + A - the matrix
2652: - lda - the leading dimension
2654: Notes:
2655: This routine is to be used in conjunction with MatSeqDenseSetPreallocation();
2656: it asserts that the preallocation has a leading dimension (the LDA parameter
2657: of Blas and Lapack fame) larger than M, the first dimension of the matrix.
2659: Level: intermediate
2661: .keywords: dense, matrix, LAPACK, BLAS
2663: .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize()
2665: @*/
2666: PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda)
2667: {
2668: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2671: 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);
2672: b->lda = lda;
2673: b->changelda = PETSC_FALSE;
2674: b->Mmax = PetscMax(b->Mmax,lda);
2675: return(0);
2676: }
2678: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2679: {
2681: PetscMPIInt size;
2684: MPI_Comm_size(comm,&size);
2685: if (size == 1) {
2686: if (scall == MAT_INITIAL_MATRIX) {
2687: MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
2688: } else {
2689: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2690: }
2691: } else {
2692: MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
2693: }
2694: return(0);
2695: }
2697: /*MC
2698: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
2700: Options Database Keys:
2701: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
2703: Level: beginner
2705: .seealso: MatCreateSeqDense()
2707: M*/
2709: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2710: {
2711: Mat_SeqDense *b;
2713: PetscMPIInt size;
2716: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
2717: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
2719: PetscNewLog(B,&b);
2720: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2721: B->data = (void*)b;
2723: b->roworiented = PETSC_TRUE;
2725: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2726: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2727: PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
2728: PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
2729: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
2730: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);
2731: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2732: #if defined(PETSC_HAVE_ELEMENTAL)
2733: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2734: #endif
2735: PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2736: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2737: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2738: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2739: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);
2740: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2741: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2742: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2743: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);
2744: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijsell_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2745: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2746: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2747: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijsell_seqdense_C",MatPtAP_SeqDense_SeqDense);
2748: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2749: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2750: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2751: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);
2753: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2754: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2755: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2756: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2757: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2758: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2759: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijsell_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2760: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2761: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2763: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2764: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2765: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2766: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
2767: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
2768: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2769: return(0);
2770: }
2772: /*@C
2773: MatDenseGetColumn - gives access to a column of a dense matrix. This is only the local part of the column. You MUST call MatDenseRestoreColumn() to avoid memory bleeding.
2775: Not Collective
2777: Input Parameter:
2778: + mat - a MATSEQDENSE or MATMPIDENSE matrix
2779: - col - column index
2781: Output Parameter:
2782: . vals - pointer to the data
2784: Level: intermediate
2786: .seealso: MatDenseRestoreColumn()
2787: @*/
2788: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
2789: {
2793: PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
2794: return(0);
2795: }
2797: /*@C
2798: MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn().
2800: Not Collective
2802: Input Parameter:
2803: . mat - a MATSEQDENSE or MATMPIDENSE matrix
2805: Output Parameter:
2806: . vals - pointer to the data
2808: Level: intermediate
2810: .seealso: MatDenseGetColumn()
2811: @*/
2812: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
2813: {
2817: PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
2818: return(0);
2819: }