Actual source code: dense.c
petsc-3.9.4 2018-09-11
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: VecGetArrayRead(x,&xx);
117: VecGetArray(b,&bb);
118: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
119: VecRestoreArrayRead(x,&xx);
120: VecRestoreArray(b,&bb);
121: }
123: for (i=0; i<N; i++) {
124: slot = l->v + rows[i]*m;
125: PetscMemzero(slot,r*sizeof(PetscScalar));
126: }
127: for (i=0; i<N; i++) {
128: slot = l->v + rows[i];
129: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
130: }
131: if (diag != 0.0) {
132: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
133: for (i=0; i<N; i++) {
134: slot = l->v + (m+1)*rows[i];
135: *slot = diag;
136: }
137: }
138: return(0);
139: }
141: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
142: {
143: Mat_SeqDense *c = (Mat_SeqDense*)(C->data);
147: MatMatMultNumeric_SeqDense_SeqDense(A,P,c->ptapwork);
148: MatTransposeMatMultNumeric_SeqDense_SeqDense(P,c->ptapwork,C);
149: return(0);
150: }
152: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat *C)
153: {
154: Mat_SeqDense *c;
158: MatCreateSeqDense(PetscObjectComm((PetscObject)A),P->cmap->N,P->cmap->N,NULL,C);
159: c = (Mat_SeqDense*)((*C)->data);
160: MatCreateSeqDense(PetscObjectComm((PetscObject)A),A->rmap->N,P->cmap->N,NULL,&c->ptapwork);
161: return(0);
162: }
164: PETSC_INTERN PetscErrorCode MatPtAP_SeqDense_SeqDense(Mat A,Mat P,MatReuse reuse,PetscReal fill,Mat *C)
165: {
169: if (reuse == MAT_INITIAL_MATRIX) {
170: MatPtAPSymbolic_SeqDense_SeqDense(A,P,fill,C);
171: }
172: PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
173: (*(*C)->ops->ptapnumeric)(A,P,*C);
174: PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
175: return(0);
176: }
178: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
179: {
180: Mat B = NULL;
181: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
182: Mat_SeqDense *b;
184: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
185: MatScalar *av=a->a;
186: PetscBool isseqdense;
189: if (reuse == MAT_REUSE_MATRIX) {
190: PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);
191: if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type);
192: }
193: if (reuse != MAT_REUSE_MATRIX) {
194: MatCreate(PetscObjectComm((PetscObject)A),&B);
195: MatSetSizes(B,m,n,m,n);
196: MatSetType(B,MATSEQDENSE);
197: MatSeqDenseSetPreallocation(B,NULL);
198: b = (Mat_SeqDense*)(B->data);
199: } else {
200: b = (Mat_SeqDense*)((*newmat)->data);
201: PetscMemzero(b->v,m*n*sizeof(PetscScalar));
202: }
203: for (i=0; i<m; i++) {
204: PetscInt j;
205: for (j=0;j<ai[1]-ai[0];j++) {
206: b->v[*aj*m+i] = *av;
207: aj++;
208: av++;
209: }
210: ai++;
211: }
213: if (reuse == MAT_INPLACE_MATRIX) {
214: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
215: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
216: MatHeaderReplace(A,&B);
217: } else {
218: if (B) *newmat = B;
219: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
220: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
221: }
222: return(0);
223: }
225: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
226: {
227: Mat B;
228: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
230: PetscInt i, j;
231: PetscInt *rows, *nnz;
232: MatScalar *aa = a->v, *vals;
235: MatCreate(PetscObjectComm((PetscObject)A),&B);
236: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
237: MatSetType(B,MATSEQAIJ);
238: PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
239: for (j=0; j<A->cmap->n; j++) {
240: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
241: aa += a->lda;
242: }
243: MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
244: aa = a->v;
245: for (j=0; j<A->cmap->n; j++) {
246: PetscInt numRows = 0;
247: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
248: MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
249: aa += a->lda;
250: }
251: PetscFree3(rows,nnz,vals);
252: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
253: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
255: if (reuse == MAT_INPLACE_MATRIX) {
256: MatHeaderReplace(A,&B);
257: } else {
258: *newmat = B;
259: }
260: return(0);
261: }
263: static PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
264: {
265: Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
266: PetscScalar oalpha = alpha;
267: PetscInt j;
268: PetscBLASInt N,m,ldax,lday,one = 1;
272: PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
273: PetscBLASIntCast(X->rmap->n,&m);
274: PetscBLASIntCast(x->lda,&ldax);
275: PetscBLASIntCast(y->lda,&lday);
276: if (ldax>m || lday>m) {
277: for (j=0; j<X->cmap->n; j++) {
278: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one));
279: }
280: } else {
281: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one));
282: }
283: PetscObjectStateIncrease((PetscObject)Y);
284: PetscLogFlops(PetscMax(2*N-1,0));
285: return(0);
286: }
288: static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
289: {
290: PetscInt N = A->rmap->n*A->cmap->n;
293: info->block_size = 1.0;
294: info->nz_allocated = (double)N;
295: info->nz_used = (double)N;
296: info->nz_unneeded = (double)0;
297: info->assemblies = (double)A->num_ass;
298: info->mallocs = 0;
299: info->memory = ((PetscObject)A)->mem;
300: info->fill_ratio_given = 0;
301: info->fill_ratio_needed = 0;
302: info->factor_mallocs = 0;
303: return(0);
304: }
306: static PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
307: {
308: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
309: PetscScalar oalpha = alpha;
311: PetscBLASInt one = 1,j,nz,lda;
314: PetscBLASIntCast(a->lda,&lda);
315: if (lda>A->rmap->n) {
316: PetscBLASIntCast(A->rmap->n,&nz);
317: for (j=0; j<A->cmap->n; j++) {
318: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one));
319: }
320: } else {
321: PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
322: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one));
323: }
324: PetscLogFlops(nz);
325: return(0);
326: }
328: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl)
329: {
330: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
331: PetscInt i,j,m = A->rmap->n,N;
332: PetscScalar *v = a->v;
335: *fl = PETSC_FALSE;
336: if (A->rmap->n != A->cmap->n) return(0);
337: N = a->lda;
339: for (i=0; i<m; i++) {
340: for (j=i+1; j<m; j++) {
341: if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) return(0);
342: }
343: }
344: *fl = PETSC_TRUE;
345: return(0);
346: }
348: static PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
349: {
350: Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l;
352: PetscInt lda = (PetscInt)mat->lda,j,m;
355: PetscLayoutReference(A->rmap,&newi->rmap);
356: PetscLayoutReference(A->cmap,&newi->cmap);
357: MatSeqDenseSetPreallocation(newi,NULL);
358: if (cpvalues == MAT_COPY_VALUES) {
359: l = (Mat_SeqDense*)newi->data;
360: if (lda>A->rmap->n) {
361: m = A->rmap->n;
362: for (j=0; j<A->cmap->n; j++) {
363: PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));
364: }
365: } else {
366: PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
367: }
368: }
369: newi->assembled = PETSC_TRUE;
370: return(0);
371: }
373: static PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
374: {
378: MatCreate(PetscObjectComm((PetscObject)A),newmat);
379: MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
380: MatSetType(*newmat,((PetscObject)A)->type_name);
381: MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
382: return(0);
383: }
386: static PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*);
388: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
389: {
390: MatFactorInfo info;
394: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
395: MatLUFactor_SeqDense(fact,0,0,&info);
396: return(0);
397: }
399: static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
400: {
401: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
402: PetscErrorCode ierr;
403: const PetscScalar *x;
404: PetscScalar *y;
405: PetscBLASInt one = 1,info,m;
408: PetscBLASIntCast(A->rmap->n,&m);
409: VecGetArrayRead(xx,&x);
410: VecGetArray(yy,&y);
411: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
412: if (A->factortype == MAT_FACTOR_LU) {
413: #if defined(PETSC_MISSING_LAPACK_GETRS)
414: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
415: #else
416: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
417: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
418: PetscFPTrapPop();
419: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
420: #endif
421: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
422: #if defined(PETSC_MISSING_LAPACK_POTRS)
423: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
424: #else
425: if (A->spd) {
426: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
427: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
428: PetscFPTrapPop();
429: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
430: #if defined(PETSC_USE_COMPLEX)
431: } else if (A->hermitian) {
432: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
433: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
434: PetscFPTrapPop();
435: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
436: #endif
437: } else { /* symmetric case */
438: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
439: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
440: PetscFPTrapPop();
441: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
442: }
443: #endif
444: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
445: VecRestoreArrayRead(xx,&x);
446: VecRestoreArray(yy,&y);
447: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
448: return(0);
449: }
451: static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
452: {
453: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
455: PetscScalar *b,*x;
456: PetscInt n;
457: PetscBLASInt nrhs,info,m;
458: PetscBool flg;
461: PetscBLASIntCast(A->rmap->n,&m);
462: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
463: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
464: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
465: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
467: MatGetSize(B,NULL,&n);
468: PetscBLASIntCast(n,&nrhs);
469: MatDenseGetArray(B,&b);
470: MatDenseGetArray(X,&x);
472: PetscMemcpy(x,b,m*nrhs*sizeof(PetscScalar));
474: if (A->factortype == MAT_FACTOR_LU) {
475: #if defined(PETSC_MISSING_LAPACK_GETRS)
476: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
477: #else
478: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
479: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
480: PetscFPTrapPop();
481: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
482: #endif
483: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
484: #if defined(PETSC_MISSING_LAPACK_POTRS)
485: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
486: #else
487: if (A->spd) {
488: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
489: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
490: PetscFPTrapPop();
491: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
492: #if defined(PETSC_USE_COMPLEX)
493: } else if (A->hermitian) {
494: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
495: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
496: PetscFPTrapPop();
497: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
498: #endif
499: } else { /* symmetric case */
500: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
501: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
502: PetscFPTrapPop();
503: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
504: }
505: #endif
506: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
508: MatDenseRestoreArray(B,&b);
509: MatDenseRestoreArray(X,&x);
510: PetscLogFlops(nrhs*(2.0*m*m - m));
511: return(0);
512: }
514: static PetscErrorCode MatConjugate_SeqDense(Mat);
516: static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
517: {
518: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
519: PetscErrorCode ierr;
520: const PetscScalar *x;
521: PetscScalar *y;
522: PetscBLASInt one = 1,info,m;
525: PetscBLASIntCast(A->rmap->n,&m);
526: VecGetArrayRead(xx,&x);
527: VecGetArray(yy,&y);
528: PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));
529: if (A->factortype == MAT_FACTOR_LU) {
530: #if defined(PETSC_MISSING_LAPACK_GETRS)
531: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable.");
532: #else
533: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
534: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
535: PetscFPTrapPop();
536: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
537: #endif
538: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
539: #if defined(PETSC_MISSING_LAPACK_POTRS)
540: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable.");
541: #else
542: if (A->spd) {
543: #if defined(PETSC_USE_COMPLEX)
544: MatConjugate_SeqDense(A);
545: #endif
546: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
547: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
548: PetscFPTrapPop();
549: #if defined(PETSC_USE_COMPLEX)
550: MatConjugate_SeqDense(A);
551: #endif
552: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
553: #if defined(PETSC_USE_COMPLEX)
554: } else if (A->hermitian) {
555: MatConjugate_SeqDense(A);
556: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
557: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
558: PetscFPTrapPop();
559: MatConjugate_SeqDense(A);
560: #endif
561: } else { /* symmetric case */
562: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
563: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
564: PetscFPTrapPop();
565: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
566: }
567: #endif
568: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
569: VecRestoreArrayRead(xx,&x);
570: VecRestoreArray(yy,&y);
571: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
572: return(0);
573: }
575: static PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
576: {
577: PetscErrorCode ierr;
578: const PetscScalar *x;
579: PetscScalar *y,sone = 1.0;
580: Vec tmp = 0;
583: VecGetArrayRead(xx,&x);
584: VecGetArray(yy,&y);
585: if (!A->rmap->n || !A->cmap->n) return(0);
586: if (yy == zz) {
587: VecDuplicate(yy,&tmp);
588: PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
589: VecCopy(yy,tmp);
590: }
591: MatSolve_SeqDense(A,xx,yy);
592: if (tmp) {
593: VecAXPY(yy,sone,tmp);
594: VecDestroy(&tmp);
595: } else {
596: VecAXPY(yy,sone,zz);
597: }
598: VecRestoreArrayRead(xx,&x);
599: VecRestoreArray(yy,&y);
600: PetscLogFlops(A->cmap->n);
601: return(0);
602: }
604: static PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
605: {
606: PetscErrorCode ierr;
607: const PetscScalar *x;
608: PetscScalar *y,sone = 1.0;
609: Vec tmp = NULL;
612: if (!A->rmap->n || !A->cmap->n) return(0);
613: VecGetArrayRead(xx,&x);
614: VecGetArray(yy,&y);
615: if (yy == zz) {
616: VecDuplicate(yy,&tmp);
617: PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);
618: VecCopy(yy,tmp);
619: }
620: MatSolveTranspose_SeqDense(A,xx,yy);
621: if (tmp) {
622: VecAXPY(yy,sone,tmp);
623: VecDestroy(&tmp);
624: } else {
625: VecAXPY(yy,sone,zz);
626: }
627: VecRestoreArrayRead(xx,&x);
628: VecRestoreArray(yy,&y);
629: PetscLogFlops(A->cmap->n);
630: return(0);
631: }
633: /* ---------------------------------------------------------------*/
634: /* COMMENT: I have chosen to hide row permutation in the pivots,
635: rather than put it in the Mat->row slot.*/
636: static PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
637: {
638: #if defined(PETSC_MISSING_LAPACK_GETRF)
640: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable.");
641: #else
642: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
644: PetscBLASInt n,m,info;
647: PetscBLASIntCast(A->cmap->n,&n);
648: PetscBLASIntCast(A->rmap->n,&m);
649: if (!mat->pivots) {
650: PetscMalloc1(A->rmap->n,&mat->pivots);
651: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
652: }
653: if (!A->rmap->n || !A->cmap->n) return(0);
654: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
655: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
656: PetscFPTrapPop();
658: if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
659: if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
661: A->ops->solve = MatSolve_SeqDense;
662: A->ops->matsolve = MatMatSolve_SeqDense;
663: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
664: A->ops->solveadd = MatSolveAdd_SeqDense;
665: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
666: A->factortype = MAT_FACTOR_LU;
668: PetscFree(A->solvertype);
669: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
671: PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
672: #endif
673: return(0);
674: }
676: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
677: static PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
678: {
679: #if defined(PETSC_MISSING_LAPACK_POTRF)
681: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
682: #else
683: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
685: PetscBLASInt info,n;
688: PetscBLASIntCast(A->cmap->n,&n);
689: if (!A->rmap->n || !A->cmap->n) return(0);
690: if (A->spd) {
691: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
692: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
693: PetscFPTrapPop();
694: #if defined(PETSC_USE_COMPLEX)
695: } else if (A->hermitian) {
696: if (!mat->pivots) {
697: PetscMalloc1(A->rmap->n,&mat->pivots);
698: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
699: }
700: if (!mat->fwork) {
701: PetscScalar dummy;
703: mat->lfwork = -1;
704: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
705: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
706: PetscFPTrapPop();
707: mat->lfwork = (PetscInt)PetscRealPart(dummy);
708: PetscMalloc1(mat->lfwork,&mat->fwork);
709: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
710: }
711: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
712: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
713: PetscFPTrapPop();
714: #endif
715: } else { /* symmetric case */
716: if (!mat->pivots) {
717: PetscMalloc1(A->rmap->n,&mat->pivots);
718: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
719: }
720: if (!mat->fwork) {
721: PetscScalar dummy;
723: mat->lfwork = -1;
724: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
725: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
726: PetscFPTrapPop();
727: mat->lfwork = (PetscInt)PetscRealPart(dummy);
728: PetscMalloc1(mat->lfwork,&mat->fwork);
729: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
730: }
731: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
732: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
733: PetscFPTrapPop();
734: }
735: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
737: A->ops->solve = MatSolve_SeqDense;
738: A->ops->matsolve = MatMatSolve_SeqDense;
739: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
740: A->ops->solveadd = MatSolveAdd_SeqDense;
741: A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
742: A->factortype = MAT_FACTOR_CHOLESKY;
744: PetscFree(A->solvertype);
745: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
747: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
748: #endif
749: return(0);
750: }
753: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
754: {
756: MatFactorInfo info;
759: info.fill = 1.0;
761: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
762: MatCholeskyFactor_SeqDense(fact,0,&info);
763: return(0);
764: }
766: static PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
767: {
769: fact->assembled = PETSC_TRUE;
770: fact->preallocated = PETSC_TRUE;
771: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
772: fact->ops->solve = MatSolve_SeqDense;
773: fact->ops->matsolve = MatMatSolve_SeqDense;
774: fact->ops->solvetranspose = MatSolveTranspose_SeqDense;
775: fact->ops->solveadd = MatSolveAdd_SeqDense;
776: fact->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
777: return(0);
778: }
780: static PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
781: {
783: fact->preallocated = PETSC_TRUE;
784: fact->assembled = PETSC_TRUE;
785: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
786: fact->ops->solve = MatSolve_SeqDense;
787: fact->ops->matsolve = MatMatSolve_SeqDense;
788: fact->ops->solvetranspose = MatSolveTranspose_SeqDense;
789: fact->ops->solveadd = MatSolveAdd_SeqDense;
790: fact->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense;
791: return(0);
792: }
794: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
795: {
799: MatCreate(PetscObjectComm((PetscObject)A),fact);
800: MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
801: MatSetType(*fact,((PetscObject)A)->type_name);
802: if (ftype == MAT_FACTOR_LU) {
803: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
804: } else {
805: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
806: }
807: (*fact)->factortype = ftype;
809: PetscFree((*fact)->solvertype);
810: PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
811: return(0);
812: }
814: /* ------------------------------------------------------------------*/
815: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
816: {
817: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
818: PetscScalar *x,*v = mat->v,zero = 0.0,xt;
819: const PetscScalar *b;
820: PetscErrorCode ierr;
821: PetscInt m = A->rmap->n,i;
822: PetscBLASInt o = 1,bm;
825: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
826: PetscBLASIntCast(m,&bm);
827: if (flag & SOR_ZERO_INITIAL_GUESS) {
828: /* this is a hack fix, should have another version without the second BLASdotu */
829: VecSet(xx,zero);
830: }
831: VecGetArray(xx,&x);
832: VecGetArrayRead(bb,&b);
833: its = its*lits;
834: if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
835: while (its--) {
836: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
837: for (i=0; i<m; i++) {
838: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
839: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
840: }
841: }
842: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
843: for (i=m-1; i>=0; i--) {
844: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
845: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
846: }
847: }
848: }
849: VecRestoreArrayRead(bb,&b);
850: VecRestoreArray(xx,&x);
851: return(0);
852: }
854: /* -----------------------------------------------------------------*/
855: PETSC_INTERN PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
856: {
857: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
858: const PetscScalar *v = mat->v,*x;
859: PetscScalar *y;
860: PetscErrorCode ierr;
861: PetscBLASInt m, n,_One=1;
862: PetscScalar _DOne=1.0,_DZero=0.0;
865: PetscBLASIntCast(A->rmap->n,&m);
866: PetscBLASIntCast(A->cmap->n,&n);
867: VecGetArrayRead(xx,&x);
868: VecGetArray(yy,&y);
869: if (!A->rmap->n || !A->cmap->n) {
870: PetscBLASInt i;
871: for (i=0; i<n; i++) y[i] = 0.0;
872: } else {
873: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
874: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
875: }
876: VecRestoreArrayRead(xx,&x);
877: VecRestoreArray(yy,&y);
878: return(0);
879: }
881: PETSC_INTERN PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
882: {
883: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
884: PetscScalar *y,_DOne=1.0,_DZero=0.0;
885: PetscErrorCode ierr;
886: PetscBLASInt m, n, _One=1;
887: const PetscScalar *v = mat->v,*x;
890: PetscBLASIntCast(A->rmap->n,&m);
891: PetscBLASIntCast(A->cmap->n,&n);
892: VecGetArrayRead(xx,&x);
893: VecGetArray(yy,&y);
894: if (!A->rmap->n || !A->cmap->n) {
895: PetscBLASInt i;
896: for (i=0; i<m; i++) y[i] = 0.0;
897: } else {
898: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
899: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
900: }
901: VecRestoreArrayRead(xx,&x);
902: VecRestoreArray(yy,&y);
903: return(0);
904: }
906: PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
907: {
908: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
909: const PetscScalar *v = mat->v,*x;
910: PetscScalar *y,_DOne=1.0;
911: PetscErrorCode ierr;
912: PetscBLASInt m, n, _One=1;
915: PetscBLASIntCast(A->rmap->n,&m);
916: PetscBLASIntCast(A->cmap->n,&n);
917: if (!A->rmap->n || !A->cmap->n) return(0);
918: if (zz != yy) {VecCopy(zz,yy);}
919: VecGetArrayRead(xx,&x);
920: VecGetArray(yy,&y);
921: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
922: VecRestoreArrayRead(xx,&x);
923: VecRestoreArray(yy,&y);
924: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
925: return(0);
926: }
928: static PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
929: {
930: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
931: const PetscScalar *v = mat->v,*x;
932: PetscScalar *y;
933: PetscErrorCode ierr;
934: PetscBLASInt m, n, _One=1;
935: PetscScalar _DOne=1.0;
938: PetscBLASIntCast(A->rmap->n,&m);
939: PetscBLASIntCast(A->cmap->n,&n);
940: if (!A->rmap->n || !A->cmap->n) return(0);
941: if (zz != yy) {VecCopy(zz,yy);}
942: VecGetArrayRead(xx,&x);
943: VecGetArray(yy,&y);
944: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
945: VecRestoreArrayRead(xx,&x);
946: VecRestoreArray(yy,&y);
947: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
948: return(0);
949: }
951: /* -----------------------------------------------------------------*/
952: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
953: {
954: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
955: PetscScalar *v;
957: PetscInt i;
960: *ncols = A->cmap->n;
961: if (cols) {
962: PetscMalloc1(A->cmap->n+1,cols);
963: for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
964: }
965: if (vals) {
966: PetscMalloc1(A->cmap->n+1,vals);
967: v = mat->v + row;
968: for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
969: }
970: return(0);
971: }
973: static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
974: {
978: if (cols) {PetscFree(*cols);}
979: if (vals) {PetscFree(*vals); }
980: return(0);
981: }
982: /* ----------------------------------------------------------------*/
983: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
984: {
985: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
986: PetscInt i,j,idx=0;
989: if (!mat->roworiented) {
990: if (addv == INSERT_VALUES) {
991: for (j=0; j<n; j++) {
992: if (indexn[j] < 0) {idx += m; continue;}
993: #if defined(PETSC_USE_DEBUG)
994: 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);
995: #endif
996: for (i=0; i<m; i++) {
997: if (indexm[i] < 0) {idx++; continue;}
998: #if defined(PETSC_USE_DEBUG)
999: 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);
1000: #endif
1001: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1002: }
1003: }
1004: } else {
1005: for (j=0; j<n; j++) {
1006: if (indexn[j] < 0) {idx += m; continue;}
1007: #if defined(PETSC_USE_DEBUG)
1008: 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);
1009: #endif
1010: for (i=0; i<m; i++) {
1011: if (indexm[i] < 0) {idx++; continue;}
1012: #if defined(PETSC_USE_DEBUG)
1013: 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);
1014: #endif
1015: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1016: }
1017: }
1018: }
1019: } else {
1020: if (addv == INSERT_VALUES) {
1021: for (i=0; i<m; i++) {
1022: if (indexm[i] < 0) { idx += n; continue;}
1023: #if defined(PETSC_USE_DEBUG)
1024: 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);
1025: #endif
1026: for (j=0; j<n; j++) {
1027: if (indexn[j] < 0) { idx++; continue;}
1028: #if defined(PETSC_USE_DEBUG)
1029: 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);
1030: #endif
1031: mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1032: }
1033: }
1034: } else {
1035: for (i=0; i<m; i++) {
1036: if (indexm[i] < 0) { idx += n; continue;}
1037: #if defined(PETSC_USE_DEBUG)
1038: 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);
1039: #endif
1040: for (j=0; j<n; j++) {
1041: if (indexn[j] < 0) { idx++; continue;}
1042: #if defined(PETSC_USE_DEBUG)
1043: 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);
1044: #endif
1045: mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1046: }
1047: }
1048: }
1049: }
1050: return(0);
1051: }
1053: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1054: {
1055: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1056: PetscInt i,j;
1059: /* row-oriented output */
1060: for (i=0; i<m; i++) {
1061: if (indexm[i] < 0) {v += n;continue;}
1062: 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);
1063: for (j=0; j<n; j++) {
1064: if (indexn[j] < 0) {v++; continue;}
1065: 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);
1066: *v++ = mat->v[indexn[j]*mat->lda + indexm[i]];
1067: }
1068: }
1069: return(0);
1070: }
1072: /* -----------------------------------------------------------------*/
1074: static PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer)
1075: {
1076: Mat_SeqDense *a;
1078: PetscInt *scols,i,j,nz,header[4];
1079: int fd;
1080: PetscMPIInt size;
1081: PetscInt *rowlengths = 0,M,N,*cols,grows,gcols;
1082: PetscScalar *vals,*svals,*v,*w;
1083: MPI_Comm comm;
1086: /* force binary viewer to load .info file if it has not yet done so */
1087: PetscViewerSetUp(viewer);
1088: PetscObjectGetComm((PetscObject)viewer,&comm);
1089: MPI_Comm_size(comm,&size);
1090: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
1091: PetscViewerBinaryGetDescriptor(viewer,&fd);
1092: PetscBinaryRead(fd,header,4,PETSC_INT);
1093: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object");
1094: M = header[1]; N = header[2]; nz = header[3];
1096: /* set global size if not set already*/
1097: if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
1098: MatSetSizes(newmat,M,N,M,N);
1099: } else {
1100: /* if sizes and type are already set, check if the vector global sizes are correct */
1101: MatGetSize(newmat,&grows,&gcols);
1102: 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);
1103: }
1104: a = (Mat_SeqDense*)newmat->data;
1105: if (!a->user_alloc) {
1106: MatSeqDenseSetPreallocation(newmat,NULL);
1107: }
1109: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */
1110: a = (Mat_SeqDense*)newmat->data;
1111: v = a->v;
1112: /* Allocate some temp space to read in the values and then flip them
1113: from row major to column major */
1114: PetscMalloc1(M*N > 0 ? M*N : 1,&w);
1115: /* read in nonzero values */
1116: PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);
1117: /* now flip the values and store them in the matrix*/
1118: for (j=0; j<N; j++) {
1119: for (i=0; i<M; i++) {
1120: *v++ =w[i*N+j];
1121: }
1122: }
1123: PetscFree(w);
1124: } else {
1125: /* read row lengths */
1126: PetscMalloc1(M+1,&rowlengths);
1127: PetscBinaryRead(fd,rowlengths,M,PETSC_INT);
1129: a = (Mat_SeqDense*)newmat->data;
1130: v = a->v;
1132: /* read column indices and nonzeros */
1133: PetscMalloc1(nz+1,&scols);
1134: cols = scols;
1135: PetscBinaryRead(fd,cols,nz,PETSC_INT);
1136: PetscMalloc1(nz+1,&svals);
1137: vals = svals;
1138: PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);
1140: /* insert into matrix */
1141: for (i=0; i<M; i++) {
1142: for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j];
1143: svals += rowlengths[i]; scols += rowlengths[i];
1144: }
1145: PetscFree(vals);
1146: PetscFree(cols);
1147: PetscFree(rowlengths);
1148: }
1149: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1150: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1151: return(0);
1152: }
1154: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1155: {
1156: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1157: PetscErrorCode ierr;
1158: PetscInt i,j;
1159: const char *name;
1160: PetscScalar *v;
1161: PetscViewerFormat format;
1162: #if defined(PETSC_USE_COMPLEX)
1163: PetscBool allreal = PETSC_TRUE;
1164: #endif
1167: PetscViewerGetFormat(viewer,&format);
1168: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1169: return(0); /* do nothing for now */
1170: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1171: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1172: for (i=0; i<A->rmap->n; i++) {
1173: v = a->v + i;
1174: PetscViewerASCIIPrintf(viewer,"row %D:",i);
1175: for (j=0; j<A->cmap->n; j++) {
1176: #if defined(PETSC_USE_COMPLEX)
1177: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1178: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1179: } else if (PetscRealPart(*v)) {
1180: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1181: }
1182: #else
1183: if (*v) {
1184: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1185: }
1186: #endif
1187: v += a->lda;
1188: }
1189: PetscViewerASCIIPrintf(viewer,"\n");
1190: }
1191: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1192: } else {
1193: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1194: #if defined(PETSC_USE_COMPLEX)
1195: /* determine if matrix has all real values */
1196: v = a->v;
1197: for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1198: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1199: }
1200: #endif
1201: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1202: PetscObjectGetName((PetscObject)A,&name);
1203: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1204: PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1205: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1206: }
1208: for (i=0; i<A->rmap->n; i++) {
1209: v = a->v + i;
1210: for (j=0; j<A->cmap->n; j++) {
1211: #if defined(PETSC_USE_COMPLEX)
1212: if (allreal) {
1213: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1214: } else {
1215: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1216: }
1217: #else
1218: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1219: #endif
1220: v += a->lda;
1221: }
1222: PetscViewerASCIIPrintf(viewer,"\n");
1223: }
1224: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1225: PetscViewerASCIIPrintf(viewer,"];\n");
1226: }
1227: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1228: }
1229: PetscViewerFlush(viewer);
1230: return(0);
1231: }
1233: static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer)
1234: {
1235: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1236: PetscErrorCode ierr;
1237: int fd;
1238: PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n;
1239: PetscScalar *v,*anonz,*vals;
1240: PetscViewerFormat format;
1243: PetscViewerBinaryGetDescriptor(viewer,&fd);
1245: PetscViewerGetFormat(viewer,&format);
1246: if (format == PETSC_VIEWER_NATIVE) {
1247: /* store the matrix as a dense matrix */
1248: PetscMalloc1(4,&col_lens);
1250: col_lens[0] = MAT_FILE_CLASSID;
1251: col_lens[1] = m;
1252: col_lens[2] = n;
1253: col_lens[3] = MATRIX_BINARY_FORMAT_DENSE;
1255: PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);
1256: PetscFree(col_lens);
1258: /* write out matrix, by rows */
1259: PetscMalloc1(m*n+1,&vals);
1260: v = a->v;
1261: for (j=0; j<n; j++) {
1262: for (i=0; i<m; i++) {
1263: vals[j + i*n] = *v++;
1264: }
1265: }
1266: PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);
1267: PetscFree(vals);
1268: } else {
1269: PetscMalloc1(4+nz,&col_lens);
1271: col_lens[0] = MAT_FILE_CLASSID;
1272: col_lens[1] = m;
1273: col_lens[2] = n;
1274: col_lens[3] = nz;
1276: /* store lengths of each row and write (including header) to file */
1277: for (i=0; i<m; i++) col_lens[4+i] = n;
1278: PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);
1280: /* Possibly should write in smaller increments, not whole matrix at once? */
1281: /* store column indices (zero start index) */
1282: ict = 0;
1283: for (i=0; i<m; i++) {
1284: for (j=0; j<n; j++) col_lens[ict++] = j;
1285: }
1286: PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);
1287: PetscFree(col_lens);
1289: /* store nonzero values */
1290: PetscMalloc1(nz+1,&anonz);
1291: ict = 0;
1292: for (i=0; i<m; i++) {
1293: v = a->v + i;
1294: for (j=0; j<n; j++) {
1295: anonz[ict++] = *v; v += a->lda;
1296: }
1297: }
1298: PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);
1299: PetscFree(anonz);
1300: }
1301: return(0);
1302: }
1304: #include <petscdraw.h>
1305: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1306: {
1307: Mat A = (Mat) Aa;
1308: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1309: PetscErrorCode ierr;
1310: PetscInt m = A->rmap->n,n = A->cmap->n,i,j;
1311: int color = PETSC_DRAW_WHITE;
1312: PetscScalar *v = a->v;
1313: PetscViewer viewer;
1314: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1315: PetscViewerFormat format;
1318: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1319: PetscViewerGetFormat(viewer,&format);
1320: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
1322: /* Loop over matrix elements drawing boxes */
1324: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1325: PetscDrawCollectiveBegin(draw);
1326: /* Blue for negative and Red for positive */
1327: for (j = 0; j < n; j++) {
1328: x_l = j; x_r = x_l + 1.0;
1329: for (i = 0; i < m; i++) {
1330: y_l = m - i - 1.0;
1331: y_r = y_l + 1.0;
1332: if (PetscRealPart(v[j*m+i]) > 0.) {
1333: color = PETSC_DRAW_RED;
1334: } else if (PetscRealPart(v[j*m+i]) < 0.) {
1335: color = PETSC_DRAW_BLUE;
1336: } else {
1337: continue;
1338: }
1339: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1340: }
1341: }
1342: PetscDrawCollectiveEnd(draw);
1343: } else {
1344: /* use contour shading to indicate magnitude of values */
1345: /* first determine max of all nonzero values */
1346: PetscReal minv = 0.0, maxv = 0.0;
1347: PetscDraw popup;
1349: for (i=0; i < m*n; i++) {
1350: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1351: }
1352: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1353: PetscDrawGetPopup(draw,&popup);
1354: PetscDrawScalePopup(popup,minv,maxv);
1356: PetscDrawCollectiveBegin(draw);
1357: for (j=0; j<n; j++) {
1358: x_l = j;
1359: x_r = x_l + 1.0;
1360: for (i=0; i<m; i++) {
1361: y_l = m - i - 1.0;
1362: y_r = y_l + 1.0;
1363: color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1364: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1365: }
1366: }
1367: PetscDrawCollectiveEnd(draw);
1368: }
1369: return(0);
1370: }
1372: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1373: {
1374: PetscDraw draw;
1375: PetscBool isnull;
1376: PetscReal xr,yr,xl,yl,h,w;
1380: PetscViewerDrawGetDraw(viewer,0,&draw);
1381: PetscDrawIsNull(draw,&isnull);
1382: if (isnull) return(0);
1384: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1385: xr += w; yr += h; xl = -w; yl = -h;
1386: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1387: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1388: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1389: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1390: PetscDrawSave(draw);
1391: return(0);
1392: }
1394: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1395: {
1397: PetscBool iascii,isbinary,isdraw;
1400: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1401: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1402: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1404: if (iascii) {
1405: MatView_SeqDense_ASCII(A,viewer);
1406: } else if (isbinary) {
1407: MatView_SeqDense_Binary(A,viewer);
1408: } else if (isdraw) {
1409: MatView_SeqDense_Draw(A,viewer);
1410: }
1411: return(0);
1412: }
1414: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar array[])
1415: {
1416: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1419: a->unplacedarray = a->v;
1420: a->unplaced_user_alloc = a->user_alloc;
1421: a->v = (PetscScalar*) array;
1422: return(0);
1423: }
1425: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1426: {
1427: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1430: a->v = a->unplacedarray;
1431: a->user_alloc = a->unplaced_user_alloc;
1432: a->unplacedarray = NULL;
1433: return(0);
1434: }
1436: static PetscErrorCode MatDestroy_SeqDense(Mat mat)
1437: {
1438: Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
1442: #if defined(PETSC_USE_LOG)
1443: PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1444: #endif
1445: PetscFree(l->pivots);
1446: PetscFree(l->fwork);
1447: MatDestroy(&l->ptapwork);
1448: if (!l->user_alloc) {PetscFree(l->v);}
1449: PetscFree(mat->data);
1451: PetscObjectChangeTypeName((PetscObject)mat,0);
1452: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1453: PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1454: PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1455: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1456: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1457: #if defined(PETSC_HAVE_ELEMENTAL)
1458: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1459: #endif
1460: PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1461: PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);
1462: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);
1463: PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);
1464: PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);
1465: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);
1466: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);
1467: PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);
1468: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
1469: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
1470: return(0);
1471: }
1473: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1474: {
1475: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1477: PetscInt k,j,m,n,M;
1478: PetscScalar *v,tmp;
1481: v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n;
1482: if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */
1483: for (j=0; j<m; j++) {
1484: for (k=0; k<j; k++) {
1485: tmp = v[j + k*M];
1486: v[j + k*M] = v[k + j*M];
1487: v[k + j*M] = tmp;
1488: }
1489: }
1490: } else { /* out-of-place transpose */
1491: Mat tmat;
1492: Mat_SeqDense *tmatd;
1493: PetscScalar *v2;
1494: PetscInt M2;
1496: if (reuse != MAT_REUSE_MATRIX) {
1497: MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1498: MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1499: MatSetType(tmat,((PetscObject)A)->type_name);
1500: MatSeqDenseSetPreallocation(tmat,NULL);
1501: } else {
1502: tmat = *matout;
1503: }
1504: tmatd = (Mat_SeqDense*)tmat->data;
1505: v = mat->v; v2 = tmatd->v; M2 = tmatd->lda;
1506: for (j=0; j<n; j++) {
1507: for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1508: }
1509: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1510: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1511: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat;
1512: else {
1513: MatHeaderMerge(A,&tmat);
1514: }
1515: }
1516: return(0);
1517: }
1519: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg)
1520: {
1521: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1522: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1523: PetscInt i,j;
1524: PetscScalar *v1,*v2;
1527: if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1528: if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1529: for (i=0; i<A1->rmap->n; i++) {
1530: v1 = mat1->v+i; v2 = mat2->v+i;
1531: for (j=0; j<A1->cmap->n; j++) {
1532: if (*v1 != *v2) {*flg = PETSC_FALSE; return(0);}
1533: v1 += mat1->lda; v2 += mat2->lda;
1534: }
1535: }
1536: *flg = PETSC_TRUE;
1537: return(0);
1538: }
1540: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1541: {
1542: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1544: PetscInt i,n,len;
1545: PetscScalar *x,zero = 0.0;
1548: VecSet(v,zero);
1549: VecGetSize(v,&n);
1550: VecGetArray(v,&x);
1551: len = PetscMin(A->rmap->n,A->cmap->n);
1552: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1553: for (i=0; i<len; i++) {
1554: x[i] = mat->v[i*mat->lda + i];
1555: }
1556: VecRestoreArray(v,&x);
1557: return(0);
1558: }
1560: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1561: {
1562: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1563: const PetscScalar *l,*r;
1564: PetscScalar x,*v;
1565: PetscErrorCode ierr;
1566: PetscInt i,j,m = A->rmap->n,n = A->cmap->n;
1569: if (ll) {
1570: VecGetSize(ll,&m);
1571: VecGetArrayRead(ll,&l);
1572: if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1573: for (i=0; i<m; i++) {
1574: x = l[i];
1575: v = mat->v + i;
1576: for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1577: }
1578: VecRestoreArrayRead(ll,&l);
1579: PetscLogFlops(1.0*n*m);
1580: }
1581: if (rr) {
1582: VecGetSize(rr,&n);
1583: VecGetArrayRead(rr,&r);
1584: if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1585: for (i=0; i<n; i++) {
1586: x = r[i];
1587: v = mat->v + i*mat->lda;
1588: for (j=0; j<m; j++) (*v++) *= x;
1589: }
1590: VecRestoreArrayRead(rr,&r);
1591: PetscLogFlops(1.0*n*m);
1592: }
1593: return(0);
1594: }
1596: static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1597: {
1598: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1599: PetscScalar *v = mat->v;
1600: PetscReal sum = 0.0;
1601: PetscInt lda =mat->lda,m=A->rmap->n,i,j;
1605: if (type == NORM_FROBENIUS) {
1606: if (lda>m) {
1607: for (j=0; j<A->cmap->n; j++) {
1608: v = mat->v+j*lda;
1609: for (i=0; i<m; i++) {
1610: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1611: }
1612: }
1613: } else {
1614: #if defined(PETSC_USE_REAL___FP16)
1615: PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1616: *nrm = BLASnrm2_(&cnt,v,&one);
1617: }
1618: #else
1619: for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1620: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1621: }
1622: }
1623: *nrm = PetscSqrtReal(sum);
1624: #endif
1625: PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1626: } else if (type == NORM_1) {
1627: *nrm = 0.0;
1628: for (j=0; j<A->cmap->n; j++) {
1629: v = mat->v + j*mat->lda;
1630: sum = 0.0;
1631: for (i=0; i<A->rmap->n; i++) {
1632: sum += PetscAbsScalar(*v); v++;
1633: }
1634: if (sum > *nrm) *nrm = sum;
1635: }
1636: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1637: } else if (type == NORM_INFINITY) {
1638: *nrm = 0.0;
1639: for (j=0; j<A->rmap->n; j++) {
1640: v = mat->v + j;
1641: sum = 0.0;
1642: for (i=0; i<A->cmap->n; i++) {
1643: sum += PetscAbsScalar(*v); v += mat->lda;
1644: }
1645: if (sum > *nrm) *nrm = sum;
1646: }
1647: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1648: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1649: return(0);
1650: }
1652: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1653: {
1654: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
1658: switch (op) {
1659: case MAT_ROW_ORIENTED:
1660: aij->roworiented = flg;
1661: break;
1662: case MAT_NEW_NONZERO_LOCATIONS:
1663: case MAT_NEW_NONZERO_LOCATION_ERR:
1664: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1665: case MAT_NEW_DIAGONALS:
1666: case MAT_KEEP_NONZERO_PATTERN:
1667: case MAT_IGNORE_OFF_PROC_ENTRIES:
1668: case MAT_USE_HASH_TABLE:
1669: case MAT_IGNORE_ZERO_ENTRIES:
1670: case MAT_IGNORE_LOWER_TRIANGULAR:
1671: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1672: break;
1673: case MAT_SPD:
1674: case MAT_SYMMETRIC:
1675: case MAT_STRUCTURALLY_SYMMETRIC:
1676: case MAT_HERMITIAN:
1677: case MAT_SYMMETRY_ETERNAL:
1678: /* These options are handled directly by MatSetOption() */
1679: break;
1680: default:
1681: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1682: }
1683: return(0);
1684: }
1686: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1687: {
1688: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1690: PetscInt lda=l->lda,m=A->rmap->n,j;
1693: if (lda>m) {
1694: for (j=0; j<A->cmap->n; j++) {
1695: PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));
1696: }
1697: } else {
1698: PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1699: }
1700: return(0);
1701: }
1703: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1704: {
1705: PetscErrorCode ierr;
1706: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1707: PetscInt m = l->lda, n = A->cmap->n, i,j;
1708: PetscScalar *slot,*bb;
1709: const PetscScalar *xx;
1712: #if defined(PETSC_USE_DEBUG)
1713: for (i=0; i<N; i++) {
1714: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1715: 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);
1716: }
1717: #endif
1719: /* fix right hand side if needed */
1720: if (x && b) {
1721: VecGetArrayRead(x,&xx);
1722: VecGetArray(b,&bb);
1723: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1724: VecRestoreArrayRead(x,&xx);
1725: VecRestoreArray(b,&bb);
1726: }
1728: for (i=0; i<N; i++) {
1729: slot = l->v + rows[i];
1730: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1731: }
1732: if (diag != 0.0) {
1733: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1734: for (i=0; i<N; i++) {
1735: slot = l->v + (m+1)*rows[i];
1736: *slot = diag;
1737: }
1738: }
1739: return(0);
1740: }
1742: static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[])
1743: {
1744: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1747: 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");
1748: *array = mat->v;
1749: return(0);
1750: }
1752: static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[])
1753: {
1755: return(0);
1756: }
1758: /*@C
1759: MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored
1761: Logically Collective on Mat
1763: Input Parameter:
1764: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1766: Output Parameter:
1767: . array - pointer to the data
1769: Level: intermediate
1771: .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
1772: @*/
1773: PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array)
1774: {
1778: PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1779: return(0);
1780: }
1782: /*@C
1783: MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
1785: Logically Collective on Mat
1787: Input Parameters:
1788: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1789: . array - pointer to the data
1791: Level: intermediate
1793: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
1794: @*/
1795: PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array)
1796: {
1800: PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1801: if (array) *array = NULL;
1802: PetscObjectStateIncrease((PetscObject)A);
1803: return(0);
1804: }
1806: /*@C
1807: MatDenseGetArrayRead - gives access to the array where the data for a SeqDense matrix is stored
1809: Not Collective
1811: Input Parameter:
1812: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1814: Output Parameter:
1815: . array - pointer to the data
1817: Level: intermediate
1819: .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead()
1820: @*/
1821: PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array)
1822: {
1826: PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
1827: return(0);
1828: }
1830: /*@C
1831: MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
1833: Not Collective
1835: Input Parameters:
1836: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1837: . array - pointer to the data
1839: Level: intermediate
1841: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray()
1842: @*/
1843: PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
1844: {
1848: PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
1849: if (array) *array = NULL;
1850: return(0);
1851: }
1853: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
1854: {
1855: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1857: PetscInt i,j,nrows,ncols;
1858: const PetscInt *irow,*icol;
1859: PetscScalar *av,*bv,*v = mat->v;
1860: Mat newmat;
1863: ISGetIndices(isrow,&irow);
1864: ISGetIndices(iscol,&icol);
1865: ISGetLocalSize(isrow,&nrows);
1866: ISGetLocalSize(iscol,&ncols);
1868: /* Check submatrixcall */
1869: if (scall == MAT_REUSE_MATRIX) {
1870: PetscInt n_cols,n_rows;
1871: MatGetSize(*B,&n_rows,&n_cols);
1872: if (n_rows != nrows || n_cols != ncols) {
1873: /* resize the result matrix to match number of requested rows/columns */
1874: MatSetSizes(*B,nrows,ncols,nrows,ncols);
1875: }
1876: newmat = *B;
1877: } else {
1878: /* Create and fill new matrix */
1879: MatCreate(PetscObjectComm((PetscObject)A),&newmat);
1880: MatSetSizes(newmat,nrows,ncols,nrows,ncols);
1881: MatSetType(newmat,((PetscObject)A)->type_name);
1882: MatSeqDenseSetPreallocation(newmat,NULL);
1883: }
1885: /* Now extract the data pointers and do the copy,column at a time */
1886: bv = ((Mat_SeqDense*)newmat->data)->v;
1888: for (i=0; i<ncols; i++) {
1889: av = v + mat->lda*icol[i];
1890: for (j=0; j<nrows; j++) *bv++ = av[irow[j]];
1891: }
1893: /* Assemble the matrices so that the correct flags are set */
1894: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
1895: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
1897: /* Free work space */
1898: ISRestoreIndices(isrow,&irow);
1899: ISRestoreIndices(iscol,&icol);
1900: *B = newmat;
1901: return(0);
1902: }
1904: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
1905: {
1907: PetscInt i;
1910: if (scall == MAT_INITIAL_MATRIX) {
1911: PetscCalloc1(n+1,B);
1912: }
1914: for (i=0; i<n; i++) {
1915: MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
1916: }
1917: return(0);
1918: }
1920: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
1921: {
1923: return(0);
1924: }
1926: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
1927: {
1929: return(0);
1930: }
1932: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
1933: {
1934: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
1936: PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;
1939: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1940: if (A->ops->copy != B->ops->copy) {
1941: MatCopy_Basic(A,B,str);
1942: return(0);
1943: }
1944: if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
1945: if (lda1>m || lda2>m) {
1946: for (j=0; j<n; j++) {
1947: PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));
1948: }
1949: } else {
1950: PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));
1951: }
1952: PetscObjectStateIncrease((PetscObject)B);
1953: return(0);
1954: }
1956: static PetscErrorCode MatSetUp_SeqDense(Mat A)
1957: {
1961: MatSeqDenseSetPreallocation(A,0);
1962: return(0);
1963: }
1965: static PetscErrorCode MatConjugate_SeqDense(Mat A)
1966: {
1967: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1968: PetscInt i,nz = A->rmap->n*A->cmap->n;
1969: PetscScalar *aa = a->v;
1972: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
1973: return(0);
1974: }
1976: static PetscErrorCode MatRealPart_SeqDense(Mat A)
1977: {
1978: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1979: PetscInt i,nz = A->rmap->n*A->cmap->n;
1980: PetscScalar *aa = a->v;
1983: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
1984: return(0);
1985: }
1987: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
1988: {
1989: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1990: PetscInt i,nz = A->rmap->n*A->cmap->n;
1991: PetscScalar *aa = a->v;
1994: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
1995: return(0);
1996: }
1998: /* ----------------------------------------------------------------*/
1999: PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2000: {
2004: if (scall == MAT_INITIAL_MATRIX) {
2005: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
2006: MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2007: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
2008: }
2009: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
2010: MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2011: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
2012: return(0);
2013: }
2015: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2016: {
2018: PetscInt m=A->rmap->n,n=B->cmap->n;
2019: Mat Cmat;
2022: 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);
2023: MatCreate(PETSC_COMM_SELF,&Cmat);
2024: MatSetSizes(Cmat,m,n,m,n);
2025: MatSetType(Cmat,MATSEQDENSE);
2026: MatSeqDenseSetPreallocation(Cmat,NULL);
2028: *C = Cmat;
2029: return(0);
2030: }
2032: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2033: {
2034: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2035: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2036: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2037: PetscBLASInt m,n,k;
2038: PetscScalar _DOne=1.0,_DZero=0.0;
2040: PetscBool flg;
2043: PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);
2044: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Second matrix must be dense");
2046: /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/
2047: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
2048: if (flg) {
2049: C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
2050: (*C->ops->matmultnumeric)(A,B,C);
2051: return(0);
2052: }
2054: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);
2055: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense");
2056: PetscBLASIntCast(C->rmap->n,&m);
2057: PetscBLASIntCast(C->cmap->n,&n);
2058: PetscBLASIntCast(A->cmap->n,&k);
2059: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2060: return(0);
2061: }
2063: PetscErrorCode MatMatTransposeMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2064: {
2068: if (scall == MAT_INITIAL_MATRIX) {
2069: PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);
2070: MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2071: PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);
2072: }
2073: PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);
2074: MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);
2075: PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);
2076: return(0);
2077: }
2079: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2080: {
2082: PetscInt m=A->rmap->n,n=B->rmap->n;
2083: Mat Cmat;
2086: 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);
2087: MatCreate(PETSC_COMM_SELF,&Cmat);
2088: MatSetSizes(Cmat,m,n,m,n);
2089: MatSetType(Cmat,MATSEQDENSE);
2090: MatSeqDenseSetPreallocation(Cmat,NULL);
2092: Cmat->assembled = PETSC_TRUE;
2094: *C = Cmat;
2095: return(0);
2096: }
2098: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2099: {
2100: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2101: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2102: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2103: PetscBLASInt m,n,k;
2104: PetscScalar _DOne=1.0,_DZero=0.0;
2108: PetscBLASIntCast(A->rmap->n,&m);
2109: PetscBLASIntCast(B->rmap->n,&n);
2110: PetscBLASIntCast(A->cmap->n,&k);
2111: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2112: return(0);
2113: }
2115: PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
2116: {
2120: if (scall == MAT_INITIAL_MATRIX) {
2121: PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);
2122: MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);
2123: PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);
2124: }
2125: PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);
2126: MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);
2127: PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);
2128: return(0);
2129: }
2131: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
2132: {
2134: PetscInt m=A->cmap->n,n=B->cmap->n;
2135: Mat Cmat;
2138: 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);
2139: MatCreate(PETSC_COMM_SELF,&Cmat);
2140: MatSetSizes(Cmat,m,n,m,n);
2141: MatSetType(Cmat,MATSEQDENSE);
2142: MatSeqDenseSetPreallocation(Cmat,NULL);
2144: Cmat->assembled = PETSC_TRUE;
2146: *C = Cmat;
2147: return(0);
2148: }
2150: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2151: {
2152: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2153: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2154: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2155: PetscBLASInt m,n,k;
2156: PetscScalar _DOne=1.0,_DZero=0.0;
2160: PetscBLASIntCast(C->rmap->n,&m);
2161: PetscBLASIntCast(C->cmap->n,&n);
2162: PetscBLASIntCast(A->rmap->n,&k);
2163: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda));
2164: return(0);
2165: }
2167: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2168: {
2169: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2171: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2172: PetscScalar *x;
2173: MatScalar *aa = a->v;
2176: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2178: VecSet(v,0.0);
2179: VecGetArray(v,&x);
2180: VecGetLocalSize(v,&p);
2181: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2182: for (i=0; i<m; i++) {
2183: x[i] = aa[i]; if (idx) idx[i] = 0;
2184: for (j=1; j<n; j++) {
2185: if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2186: }
2187: }
2188: VecRestoreArray(v,&x);
2189: return(0);
2190: }
2192: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2193: {
2194: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2196: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2197: PetscScalar *x;
2198: PetscReal atmp;
2199: MatScalar *aa = a->v;
2202: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2204: VecSet(v,0.0);
2205: VecGetArray(v,&x);
2206: VecGetLocalSize(v,&p);
2207: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2208: for (i=0; i<m; i++) {
2209: x[i] = PetscAbsScalar(aa[i]);
2210: for (j=1; j<n; j++) {
2211: atmp = PetscAbsScalar(aa[i+m*j]);
2212: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2213: }
2214: }
2215: VecRestoreArray(v,&x);
2216: return(0);
2217: }
2219: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2220: {
2221: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2223: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2224: PetscScalar *x;
2225: MatScalar *aa = a->v;
2228: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2230: VecSet(v,0.0);
2231: VecGetArray(v,&x);
2232: VecGetLocalSize(v,&p);
2233: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2234: for (i=0; i<m; i++) {
2235: x[i] = aa[i]; if (idx) idx[i] = 0;
2236: for (j=1; j<n; j++) {
2237: if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;}
2238: }
2239: }
2240: VecRestoreArray(v,&x);
2241: return(0);
2242: }
2244: static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2245: {
2246: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2248: PetscScalar *x;
2251: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2253: VecGetArray(v,&x);
2254: PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));
2255: VecRestoreArray(v,&x);
2256: return(0);
2257: }
2259: PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2260: {
2262: PetscInt i,j,m,n;
2263: PetscScalar *a;
2266: MatGetSize(A,&m,&n);
2267: PetscMemzero(norms,n*sizeof(PetscReal));
2268: MatDenseGetArray(A,&a);
2269: if (type == NORM_2) {
2270: for (i=0; i<n; i++) {
2271: for (j=0; j<m; j++) {
2272: norms[i] += PetscAbsScalar(a[j]*a[j]);
2273: }
2274: a += m;
2275: }
2276: } else if (type == NORM_1) {
2277: for (i=0; i<n; i++) {
2278: for (j=0; j<m; j++) {
2279: norms[i] += PetscAbsScalar(a[j]);
2280: }
2281: a += m;
2282: }
2283: } else if (type == NORM_INFINITY) {
2284: for (i=0; i<n; i++) {
2285: for (j=0; j<m; j++) {
2286: norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2287: }
2288: a += m;
2289: }
2290: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2291: MatDenseRestoreArray(A,&a);
2292: if (type == NORM_2) {
2293: for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2294: }
2295: return(0);
2296: }
2298: static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2299: {
2301: PetscScalar *a;
2302: PetscInt m,n,i;
2305: MatGetSize(x,&m,&n);
2306: MatDenseGetArray(x,&a);
2307: for (i=0; i<m*n; i++) {
2308: PetscRandomGetValue(rctx,a+i);
2309: }
2310: MatDenseRestoreArray(x,&a);
2311: return(0);
2312: }
2314: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d)
2315: {
2317: *missing = PETSC_FALSE;
2318: return(0);
2319: }
2321: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2322: {
2323: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2326: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2327: *vals = a->v+col*a->lda;
2328: return(0);
2329: }
2331: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2332: {
2334: *vals = 0; /* user cannot accidently use the array later */
2335: return(0);
2336: }
2338: /* -------------------------------------------------------------------*/
2339: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2340: MatGetRow_SeqDense,
2341: MatRestoreRow_SeqDense,
2342: MatMult_SeqDense,
2343: /* 4*/ MatMultAdd_SeqDense,
2344: MatMultTranspose_SeqDense,
2345: MatMultTransposeAdd_SeqDense,
2346: 0,
2347: 0,
2348: 0,
2349: /* 10*/ 0,
2350: MatLUFactor_SeqDense,
2351: MatCholeskyFactor_SeqDense,
2352: MatSOR_SeqDense,
2353: MatTranspose_SeqDense,
2354: /* 15*/ MatGetInfo_SeqDense,
2355: MatEqual_SeqDense,
2356: MatGetDiagonal_SeqDense,
2357: MatDiagonalScale_SeqDense,
2358: MatNorm_SeqDense,
2359: /* 20*/ MatAssemblyBegin_SeqDense,
2360: MatAssemblyEnd_SeqDense,
2361: MatSetOption_SeqDense,
2362: MatZeroEntries_SeqDense,
2363: /* 24*/ MatZeroRows_SeqDense,
2364: 0,
2365: 0,
2366: 0,
2367: 0,
2368: /* 29*/ MatSetUp_SeqDense,
2369: 0,
2370: 0,
2371: 0,
2372: 0,
2373: /* 34*/ MatDuplicate_SeqDense,
2374: 0,
2375: 0,
2376: 0,
2377: 0,
2378: /* 39*/ MatAXPY_SeqDense,
2379: MatCreateSubMatrices_SeqDense,
2380: 0,
2381: MatGetValues_SeqDense,
2382: MatCopy_SeqDense,
2383: /* 44*/ MatGetRowMax_SeqDense,
2384: MatScale_SeqDense,
2385: MatShift_Basic,
2386: 0,
2387: MatZeroRowsColumns_SeqDense,
2388: /* 49*/ MatSetRandom_SeqDense,
2389: 0,
2390: 0,
2391: 0,
2392: 0,
2393: /* 54*/ 0,
2394: 0,
2395: 0,
2396: 0,
2397: 0,
2398: /* 59*/ 0,
2399: MatDestroy_SeqDense,
2400: MatView_SeqDense,
2401: 0,
2402: 0,
2403: /* 64*/ 0,
2404: 0,
2405: 0,
2406: 0,
2407: 0,
2408: /* 69*/ MatGetRowMaxAbs_SeqDense,
2409: 0,
2410: 0,
2411: 0,
2412: 0,
2413: /* 74*/ 0,
2414: 0,
2415: 0,
2416: 0,
2417: 0,
2418: /* 79*/ 0,
2419: 0,
2420: 0,
2421: 0,
2422: /* 83*/ MatLoad_SeqDense,
2423: 0,
2424: MatIsHermitian_SeqDense,
2425: 0,
2426: 0,
2427: 0,
2428: /* 89*/ MatMatMult_SeqDense_SeqDense,
2429: MatMatMultSymbolic_SeqDense_SeqDense,
2430: MatMatMultNumeric_SeqDense_SeqDense,
2431: MatPtAP_SeqDense_SeqDense,
2432: MatPtAPSymbolic_SeqDense_SeqDense,
2433: /* 94*/ MatPtAPNumeric_SeqDense_SeqDense,
2434: MatMatTransposeMult_SeqDense_SeqDense,
2435: MatMatTransposeMultSymbolic_SeqDense_SeqDense,
2436: MatMatTransposeMultNumeric_SeqDense_SeqDense,
2437: 0,
2438: /* 99*/ 0,
2439: 0,
2440: 0,
2441: MatConjugate_SeqDense,
2442: 0,
2443: /*104*/ 0,
2444: MatRealPart_SeqDense,
2445: MatImaginaryPart_SeqDense,
2446: 0,
2447: 0,
2448: /*109*/ 0,
2449: 0,
2450: MatGetRowMin_SeqDense,
2451: MatGetColumnVector_SeqDense,
2452: MatMissingDiagonal_SeqDense,
2453: /*114*/ 0,
2454: 0,
2455: 0,
2456: 0,
2457: 0,
2458: /*119*/ 0,
2459: 0,
2460: 0,
2461: 0,
2462: 0,
2463: /*124*/ 0,
2464: MatGetColumnNorms_SeqDense,
2465: 0,
2466: 0,
2467: 0,
2468: /*129*/ 0,
2469: MatTransposeMatMult_SeqDense_SeqDense,
2470: MatTransposeMatMultSymbolic_SeqDense_SeqDense,
2471: MatTransposeMatMultNumeric_SeqDense_SeqDense,
2472: 0,
2473: /*134*/ 0,
2474: 0,
2475: 0,
2476: 0,
2477: 0,
2478: /*139*/ 0,
2479: 0,
2480: 0,
2481: 0,
2482: 0,
2483: /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqDense
2484: };
2486: /*@C
2487: MatCreateSeqDense - Creates a sequential dense matrix that
2488: is stored in column major order (the usual Fortran 77 manner). Many
2489: of the matrix operations use the BLAS and LAPACK routines.
2491: Collective on MPI_Comm
2493: Input Parameters:
2494: + comm - MPI communicator, set to PETSC_COMM_SELF
2495: . m - number of rows
2496: . n - number of columns
2497: - data - optional location of matrix data in column major order. Set data=NULL for PETSc
2498: to control all matrix memory allocation.
2500: Output Parameter:
2501: . A - the matrix
2503: Notes:
2504: The data input variable is intended primarily for Fortran programmers
2505: who wish to allocate their own matrix memory space. Most users should
2506: set data=NULL.
2508: Level: intermediate
2510: .keywords: dense, matrix, LAPACK, BLAS
2512: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2513: @*/
2514: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2515: {
2519: MatCreate(comm,A);
2520: MatSetSizes(*A,m,n,m,n);
2521: MatSetType(*A,MATSEQDENSE);
2522: MatSeqDenseSetPreallocation(*A,data);
2523: return(0);
2524: }
2526: /*@C
2527: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
2529: Collective on MPI_Comm
2531: Input Parameters:
2532: + B - the matrix
2533: - data - the array (or NULL)
2535: Notes:
2536: The data input variable is intended primarily for Fortran programmers
2537: who wish to allocate their own matrix memory space. Most users should
2538: need not call this routine.
2540: Level: intermediate
2542: .keywords: dense, matrix, LAPACK, BLAS
2544: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA()
2546: @*/
2547: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2548: {
2552: PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2553: return(0);
2554: }
2556: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2557: {
2558: Mat_SeqDense *b;
2562: B->preallocated = PETSC_TRUE;
2564: PetscLayoutSetUp(B->rmap);
2565: PetscLayoutSetUp(B->cmap);
2567: b = (Mat_SeqDense*)B->data;
2568: b->Mmax = B->rmap->n;
2569: b->Nmax = B->cmap->n;
2570: if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n;
2572: PetscIntMultError(b->lda,b->Nmax,NULL);
2573: if (!data) { /* petsc-allocated storage */
2574: if (!b->user_alloc) { PetscFree(b->v); }
2575: PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);
2576: PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));
2578: b->user_alloc = PETSC_FALSE;
2579: } else { /* user-allocated storage */
2580: if (!b->user_alloc) { PetscFree(b->v); }
2581: b->v = data;
2582: b->user_alloc = PETSC_TRUE;
2583: }
2584: B->assembled = PETSC_TRUE;
2585: return(0);
2586: }
2588: #if defined(PETSC_HAVE_ELEMENTAL)
2589: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2590: {
2591: Mat mat_elemental;
2593: PetscScalar *array,*v_colwise;
2594: PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;
2597: PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2598: MatDenseGetArray(A,&array);
2599: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2600: k = 0;
2601: for (j=0; j<N; j++) {
2602: cols[j] = j;
2603: for (i=0; i<M; i++) {
2604: v_colwise[j*M+i] = array[k++];
2605: }
2606: }
2607: for (i=0; i<M; i++) {
2608: rows[i] = i;
2609: }
2610: MatDenseRestoreArray(A,&array);
2612: MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2613: MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2614: MatSetType(mat_elemental,MATELEMENTAL);
2615: MatSetUp(mat_elemental);
2617: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2618: MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2619: MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2620: MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2621: PetscFree3(v_colwise,rows,cols);
2623: if (reuse == MAT_INPLACE_MATRIX) {
2624: MatHeaderReplace(A,&mat_elemental);
2625: } else {
2626: *newmat = mat_elemental;
2627: }
2628: return(0);
2629: }
2630: #endif
2632: /*@C
2633: MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array
2635: Input parameter:
2636: + A - the matrix
2637: - lda - the leading dimension
2639: Notes:
2640: This routine is to be used in conjunction with MatSeqDenseSetPreallocation();
2641: it asserts that the preallocation has a leading dimension (the LDA parameter
2642: of Blas and Lapack fame) larger than M, the first dimension of the matrix.
2644: Level: intermediate
2646: .keywords: dense, matrix, LAPACK, BLAS
2648: .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize()
2650: @*/
2651: PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda)
2652: {
2653: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2656: 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);
2657: b->lda = lda;
2658: b->changelda = PETSC_FALSE;
2659: b->Mmax = PetscMax(b->Mmax,lda);
2660: return(0);
2661: }
2663: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2664: {
2666: PetscMPIInt size;
2669: MPI_Comm_size(comm,&size);
2670: if (size == 1) {
2671: if (scall == MAT_INITIAL_MATRIX) {
2672: MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
2673: } else {
2674: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2675: }
2676: } else {
2677: MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
2678: }
2679: return(0);
2680: }
2682: /*MC
2683: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
2685: Options Database Keys:
2686: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
2688: Level: beginner
2690: .seealso: MatCreateSeqDense()
2692: M*/
2694: PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B)
2695: {
2696: Mat_SeqDense *b;
2698: PetscMPIInt size;
2701: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
2702: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
2704: PetscNewLog(B,&b);
2705: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
2706: B->data = (void*)b;
2708: b->roworiented = PETSC_TRUE;
2710: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
2711: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
2712: PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
2713: PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
2714: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
2715: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);
2716: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
2717: #if defined(PETSC_HAVE_ELEMENTAL)
2718: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
2719: #endif
2720: PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
2721: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2722: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2723: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2724: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);
2725: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2726: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2727: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2728: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);
2729: PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);
2730: PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);
2731: PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);
2732: PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);
2734: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2735: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2736: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2737: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2738: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2739: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2741: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);
2742: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);
2743: PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);
2744: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
2745: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
2746: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
2747: return(0);
2748: }
2750: /*@C
2751: 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.
2753: Not Collective
2755: Input Parameter:
2756: + mat - a MATSEQDENSE or MATMPIDENSE matrix
2757: - col - column index
2759: Output Parameter:
2760: . vals - pointer to the data
2762: Level: intermediate
2764: .seealso: MatDenseRestoreColumn()
2765: @*/
2766: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
2767: {
2771: PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
2772: return(0);
2773: }
2775: /*@C
2776: MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn().
2778: Not Collective
2780: Input Parameter:
2781: . mat - a MATSEQDENSE or MATMPIDENSE matrix
2783: Output Parameter:
2784: . vals - pointer to the data
2786: Level: intermediate
2788: .seealso: MatDenseGetColumn()
2789: @*/
2790: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
2791: {
2795: PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
2796: return(0);
2797: }