Actual source code: dense.c
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: PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian)
12: {
13: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
14: PetscInt j, k, n = A->rmap->n;
15: PetscScalar *v;
19: if (A->rmap->n != A->cmap->n) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot symmetrize a rectangular matrix");
20: MatDenseGetArray(A,&v);
21: if (!hermitian) {
22: for (k=0;k<n;k++) {
23: for (j=k;j<n;j++) {
24: v[j*mat->lda + k] = v[k*mat->lda + j];
25: }
26: }
27: } else {
28: for (k=0;k<n;k++) {
29: for (j=k;j<n;j++) {
30: v[j*mat->lda + k] = PetscConj(v[k*mat->lda + j]);
31: }
32: }
33: }
34: MatDenseRestoreArray(A,&v);
35: return(0);
36: }
38: PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A)
39: {
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);
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);
83: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
85: A->ops->solve = NULL;
86: A->ops->matsolve = NULL;
87: A->ops->solvetranspose = NULL;
88: A->ops->matsolvetranspose = NULL;
89: A->ops->solveadd = NULL;
90: A->ops->solvetransposeadd = NULL;
91: A->factortype = MAT_FACTOR_NONE;
92: PetscFree(A->solvertype);
93: return(0);
94: }
96: PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
97: {
98: PetscErrorCode ierr;
99: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
100: PetscInt m = l->lda, n = A->cmap->n,r = A->rmap->n, i,j;
101: PetscScalar *slot,*bb,*v;
102: const PetscScalar *xx;
105: if (PetscDefined(USE_DEBUG)) {
106: for (i=0; i<N; i++) {
107: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
108: 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);
109: 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);
110: }
111: }
112: if (!N) return(0);
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: MatDenseGetArray(A,&v);
132: for (i=0; i<N; i++) {
133: slot = v + rows[i]*m;
134: PetscArrayzero(slot,r);
135: }
136: for (i=0; i<N; i++) {
137: slot = v + rows[i];
138: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
139: }
140: if (diag != 0.0) {
141: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
142: for (i=0; i<N; i++) {
143: slot = v + (m+1)*rows[i];
144: *slot = diag;
145: }
146: }
147: MatDenseRestoreArray(A,&v);
148: return(0);
149: }
151: PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
152: {
153: Mat_SeqDense *c = (Mat_SeqDense*)(C->data);
157: if (c->ptapwork) {
158: (*C->ops->matmultnumeric)(A,P,c->ptapwork);
159: (*C->ops->transposematmultnumeric)(P,c->ptapwork,C);
160: } else SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Must call MatPtAPSymbolic_SeqDense_SeqDense() first");
161: return(0);
162: }
164: PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat C)
165: {
166: Mat_SeqDense *c;
167: PetscBool cisdense;
171: MatSetSizes(C,P->cmap->n,P->cmap->n,P->cmap->N,P->cmap->N);
172: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
173: if (!cisdense) {
174: PetscBool flg;
176: PetscObjectTypeCompare((PetscObject)P,((PetscObject)A)->type_name,&flg);
177: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
178: }
179: MatSetUp(C);
180: c = (Mat_SeqDense*)C->data;
181: MatCreate(PetscObjectComm((PetscObject)A),&c->ptapwork);
182: MatSetSizes(c->ptapwork,A->rmap->n,P->cmap->n,A->rmap->N,P->cmap->N);
183: MatSetType(c->ptapwork,((PetscObject)C)->type_name);
184: MatSetUp(c->ptapwork);
185: return(0);
186: }
188: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
189: {
190: Mat B = NULL;
191: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
192: Mat_SeqDense *b;
193: PetscErrorCode ierr;
194: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
195: const MatScalar *av;
196: PetscBool isseqdense;
199: if (reuse == MAT_REUSE_MATRIX) {
200: PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);
201: if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type_name);
202: }
203: if (reuse != MAT_REUSE_MATRIX) {
204: MatCreate(PetscObjectComm((PetscObject)A),&B);
205: MatSetSizes(B,m,n,m,n);
206: MatSetType(B,MATSEQDENSE);
207: MatSeqDenseSetPreallocation(B,NULL);
208: b = (Mat_SeqDense*)(B->data);
209: } else {
210: b = (Mat_SeqDense*)((*newmat)->data);
211: PetscArrayzero(b->v,m*n);
212: }
213: MatSeqAIJGetArrayRead(A,&av);
214: for (i=0; i<m; i++) {
215: PetscInt j;
216: for (j=0;j<ai[1]-ai[0];j++) {
217: b->v[*aj*m+i] = *av;
218: aj++;
219: av++;
220: }
221: ai++;
222: }
223: MatSeqAIJRestoreArrayRead(A,&av);
225: if (reuse == MAT_INPLACE_MATRIX) {
226: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
227: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
228: MatHeaderReplace(A,&B);
229: } else {
230: if (B) *newmat = B;
231: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
232: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
233: }
234: return(0);
235: }
237: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
238: {
239: Mat B = NULL;
240: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
242: PetscInt i, j;
243: PetscInt *rows, *nnz;
244: MatScalar *aa = a->v, *vals;
247: PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
248: if (reuse != MAT_REUSE_MATRIX) {
249: MatCreate(PetscObjectComm((PetscObject)A),&B);
250: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
251: MatSetType(B,MATSEQAIJ);
252: for (j=0; j<A->cmap->n; j++) {
253: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) ++nnz[i];
254: aa += a->lda;
255: }
256: MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
257: } else B = *newmat;
258: aa = a->v;
259: for (j=0; j<A->cmap->n; j++) {
260: PetscInt numRows = 0;
261: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) {rows[numRows] = i; vals[numRows++] = aa[i];}
262: MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
263: aa += a->lda;
264: }
265: PetscFree3(rows,nnz,vals);
266: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
267: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
269: if (reuse == MAT_INPLACE_MATRIX) {
270: MatHeaderReplace(A,&B);
271: } else if (reuse != MAT_REUSE_MATRIX) *newmat = B;
272: return(0);
273: }
275: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
276: {
277: Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
278: const PetscScalar *xv;
279: PetscScalar *yv;
280: PetscBLASInt N,m,ldax = 0,lday = 0,one = 1;
281: PetscErrorCode ierr;
284: MatDenseGetArrayRead(X,&xv);
285: MatDenseGetArray(Y,&yv);
286: PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
287: PetscBLASIntCast(X->rmap->n,&m);
288: PetscBLASIntCast(x->lda,&ldax);
289: PetscBLASIntCast(y->lda,&lday);
290: if (ldax>m || lday>m) {
291: PetscInt j;
293: for (j=0; j<X->cmap->n; j++) {
294: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,xv+j*ldax,&one,yv+j*lday,&one));
295: }
296: } else {
297: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&alpha,xv,&one,yv,&one));
298: }
299: MatDenseRestoreArrayRead(X,&xv);
300: MatDenseRestoreArray(Y,&yv);
301: PetscLogFlops(PetscMax(2.0*N-1,0));
302: return(0);
303: }
305: static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
306: {
307: PetscLogDouble N = A->rmap->n*A->cmap->n;
310: info->block_size = 1.0;
311: info->nz_allocated = N;
312: info->nz_used = N;
313: info->nz_unneeded = 0;
314: info->assemblies = A->num_ass;
315: info->mallocs = 0;
316: info->memory = ((PetscObject)A)->mem;
317: info->fill_ratio_given = 0;
318: info->fill_ratio_needed = 0;
319: info->factor_mallocs = 0;
320: return(0);
321: }
323: PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
324: {
325: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
326: PetscScalar *v;
328: PetscBLASInt one = 1,j,nz,lda = 0;
331: MatDenseGetArray(A,&v);
332: PetscBLASIntCast(a->lda,&lda);
333: if (lda>A->rmap->n) {
334: PetscBLASIntCast(A->rmap->n,&nz);
335: for (j=0; j<A->cmap->n; j++) {
336: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v+j*lda,&one));
337: }
338: } else {
339: PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
340: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v,&one));
341: }
342: PetscLogFlops(nz);
343: MatDenseRestoreArray(A,&v);
344: return(0);
345: }
347: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl)
348: {
349: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
350: PetscInt i,j,m = A->rmap->n,N = a->lda;
351: const PetscScalar *v;
352: PetscErrorCode ierr;
355: *fl = PETSC_FALSE;
356: if (A->rmap->n != A->cmap->n) return(0);
357: MatDenseGetArrayRead(A,&v);
358: for (i=0; i<m; i++) {
359: for (j=i; j<m; j++) {
360: if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) {
361: goto restore;
362: }
363: }
364: }
365: *fl = PETSC_TRUE;
366: restore:
367: MatDenseRestoreArrayRead(A,&v);
368: return(0);
369: }
371: static PetscErrorCode MatIsSymmetric_SeqDense(Mat A,PetscReal rtol,PetscBool *fl)
372: {
373: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
374: PetscInt i,j,m = A->rmap->n,N = a->lda;
375: const PetscScalar *v;
376: PetscErrorCode ierr;
379: *fl = PETSC_FALSE;
380: if (A->rmap->n != A->cmap->n) return(0);
381: MatDenseGetArrayRead(A,&v);
382: for (i=0; i<m; i++) {
383: for (j=i; j<m; j++) {
384: if (PetscAbsScalar(v[i+j*N] - v[j+i*N]) > rtol) {
385: goto restore;
386: }
387: }
388: }
389: *fl = PETSC_TRUE;
390: restore:
391: MatDenseRestoreArrayRead(A,&v);
392: return(0);
393: }
395: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
396: {
397: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
399: PetscInt lda = (PetscInt)mat->lda,j,m,nlda = lda;
402: PetscLayoutReference(A->rmap,&newi->rmap);
403: PetscLayoutReference(A->cmap,&newi->cmap);
404: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */
405: MatDenseSetLDA(newi,lda);
406: }
407: MatSeqDenseSetPreallocation(newi,NULL);
408: if (cpvalues == MAT_COPY_VALUES) {
409: const PetscScalar *av;
410: PetscScalar *v;
412: MatDenseGetArrayRead(A,&av);
413: MatDenseGetArray(newi,&v);
414: MatDenseGetLDA(newi,&nlda);
415: m = A->rmap->n;
416: if (lda>m || nlda>m) {
417: for (j=0; j<A->cmap->n; j++) {
418: PetscArraycpy(v+j*nlda,av+j*lda,m);
419: }
420: } else {
421: PetscArraycpy(v,av,A->rmap->n*A->cmap->n);
422: }
423: MatDenseRestoreArray(newi,&v);
424: MatDenseRestoreArrayRead(A,&av);
425: }
426: return(0);
427: }
429: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
430: {
434: MatCreate(PetscObjectComm((PetscObject)A),newmat);
435: MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
436: MatSetType(*newmat,((PetscObject)A)->type_name);
437: MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
438: return(0);
439: }
441: static PetscErrorCode MatSolve_SeqDense_Internal_LU(Mat A, PetscScalar *x, PetscBLASInt xlda, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
442: {
443: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
444: PetscBLASInt info;
445: PetscErrorCode ierr;
448: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
449: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_(T ? "T" : "N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
450: PetscFPTrapPop();
451: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
452: PetscLogFlops(nrhs*(2.0*m*m - m));
453: return(0);
454: }
456: static PetscErrorCode MatConjugate_SeqDense(Mat);
458: static PetscErrorCode MatSolve_SeqDense_Internal_Cholesky(Mat A, PetscScalar *x, PetscBLASInt xlda, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
459: {
460: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
461: PetscBLASInt info;
462: PetscErrorCode ierr;
465: if (A->spd) {
466: if (PetscDefined(USE_COMPLEX) && T) {MatConjugate_SeqDense(A);}
467: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
468: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
469: PetscFPTrapPop();
470: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
471: if (PetscDefined(USE_COMPLEX) && T) {MatConjugate_SeqDense(A);}
472: #if defined(PETSC_USE_COMPLEX)
473: } else if (A->hermitian) {
474: if (T) {MatConjugate_SeqDense(A);}
475: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
476: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
477: PetscFPTrapPop();
478: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
479: if (T) {MatConjugate_SeqDense(A);}
480: #endif
481: } else { /* symmetric case */
482: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
483: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
484: PetscFPTrapPop();
485: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
486: }
487: PetscLogFlops(nrhs*(2.0*m*m - m));
488: return(0);
489: }
491: static PetscErrorCode MatSolve_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt xlda, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
492: {
493: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
494: PetscBLASInt info;
495: char trans;
496: PetscErrorCode ierr;
499: if (PetscDefined(USE_COMPLEX)) {
500: trans = 'C';
501: } else {
502: trans = 'T';
503: }
504: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
505: PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("L", &trans, &m,&nrhs,&mat->rank,mat->v,&mat->lda,mat->tau,x,&xlda,mat->fwork,&mat->lfwork,&info));
506: PetscFPTrapPop();
507: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ORMQR - Bad orthogonal transform");
508: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
509: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U", "N", "N", &mat->rank,&nrhs,mat->v,&mat->lda,x,&xlda,&info));
510: PetscFPTrapPop();
511: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"TRTRS - Bad triangular solve");
512: for (PetscInt j = 0; j < nrhs; j++) {
513: for (PetscInt i = mat->rank; i < k; i++) {
514: x[j*xlda + i] = 0.;
515: }
516: }
517: PetscLogFlops(nrhs*(4.0*m*mat->rank - PetscSqr(mat->rank)));
518: return(0);
519: }
521: static PetscErrorCode MatSolveTranspose_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt xlda, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
522: {
523: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
524: PetscBLASInt info;
525: PetscErrorCode ierr;
528: if (A->rmap->n == A->cmap->n && mat->rank == A->rmap->n) {
529: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
530: PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U", "T", "N", &m,&nrhs,mat->v,&mat->lda,x,&xlda,&info));
531: PetscFPTrapPop();
532: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"TRTRS - Bad triangular solve");
533: if (PetscDefined(USE_COMPLEX)) {MatConjugate_SeqDense(A);}
534: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
535: PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("L", "N", &m,&nrhs,&mat->rank,mat->v,&mat->lda,mat->tau,x,&xlda,mat->fwork,&mat->lfwork,&info));
536: PetscFPTrapPop();
537: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ORMQR - Bad orthogonal transform");
538: if (PetscDefined(USE_COMPLEX)) {MatConjugate_SeqDense(A);}
539: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"QR factored matrix cannot be used for transpose solve");
540: PetscLogFlops(nrhs*(4.0*m*mat->rank - PetscSqr(mat->rank)));
541: return(0);
542: }
544: static PetscErrorCode MatSolve_SeqDense_SetUp(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
545: {
546: Mat_SeqDense *mat = (Mat_SeqDense *) A->data;
547: PetscScalar *y;
548: PetscBLASInt m=0, k=0;
549: PetscErrorCode ierr;
552: PetscBLASIntCast(A->rmap->n,&m);
553: PetscBLASIntCast(A->cmap->n,&k);
554: if (k < m) {
555: VecCopy(xx, mat->qrrhs);
556: VecGetArray(mat->qrrhs,&y);
557: } else {
558: VecCopy(xx, yy);
559: VecGetArray(yy,&y);
560: }
561: *_y = y;
562: *_k = k;
563: *_m = m;
564: return(0);
565: }
567: static PetscErrorCode MatSolve_SeqDense_TearDown(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
568: {
569: Mat_SeqDense *mat = (Mat_SeqDense *) A->data;
570: PetscScalar *y = NULL;
571: PetscBLASInt m, k;
575: y = *_y;
576: *_y = NULL;
577: k = *_k;
578: m = *_m;
579: if (k < m) {
580: PetscScalar *yv;
581: VecGetArray(yy,&yv);
582: PetscArraycpy(yv, y, k);
583: VecRestoreArray(yy,&yv);
584: VecRestoreArray(mat->qrrhs, &y);
585: } else {
586: VecRestoreArray(yy,&y);
587: }
588: return(0);
589: }
591: static PetscErrorCode MatSolve_SeqDense_LU(Mat A, Vec xx, Vec yy)
592: {
593: PetscScalar *y = NULL;
594: PetscBLASInt m = 0, k = 0;
598: MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
599: MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_FALSE);
600: MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
601: return(0);
602: }
604: static PetscErrorCode MatSolveTranspose_SeqDense_LU(Mat A, Vec xx, Vec yy)
605: {
606: PetscScalar *y = NULL;
607: PetscBLASInt m = 0, k = 0;
611: MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
612: MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_TRUE);
613: MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
614: return(0);
615: }
617: static PetscErrorCode MatSolve_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
618: {
619: PetscScalar *y = NULL;
620: PetscBLASInt m = 0, k = 0;
624: MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
625: MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_FALSE);
626: MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
627: return(0);
628: }
630: static PetscErrorCode MatSolveTranspose_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
631: {
632: PetscScalar *y = NULL;
633: PetscBLASInt m = 0, k = 0;
637: MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
638: MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_TRUE);
639: MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
640: return(0);
641: }
643: static PetscErrorCode MatSolve_SeqDense_QR(Mat A, Vec xx, Vec yy)
644: {
645: PetscScalar *y = NULL;
646: PetscBLASInt m = 0, k = 0;
650: MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
651: MatSolve_SeqDense_Internal_QR(A, y, PetscMax(m,k), m, 1, k);
652: MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
653: return(0);
654: }
656: static PetscErrorCode MatSolveTranspose_SeqDense_QR(Mat A, Vec xx, Vec yy)
657: {
658: PetscScalar *y = NULL;
659: PetscBLASInt m = 0, k = 0;
663: MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);
664: MatSolveTranspose_SeqDense_Internal_QR(A, y, PetscMax(m,k), m, 1, k);
665: MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);
666: return(0);
667: }
669: static PetscErrorCode MatMatSolve_SeqDense_SetUp(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ylda, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
670: {
671: PetscErrorCode ierr;
672: const PetscScalar *b;
673: PetscScalar *y;
674: PetscInt n, _blda, _xlda;
675: PetscBLASInt nrhs=0,m=0,k=0,blda=0,xlda=0,ylda=0;
678: *_ylda=0; *_m=0; *_nrhs=0; *_k=0;
679: PetscBLASIntCast(A->rmap->n,&m);
680: PetscBLASIntCast(A->cmap->n,&k);
681: MatGetSize(B,NULL,&n);
682: PetscBLASIntCast(n,&nrhs);
683: MatDenseGetLDA(B,&_blda);
684: PetscBLASIntCast(_blda, &blda);
685: MatDenseGetLDA(X,&_xlda);
686: PetscBLASIntCast(_xlda, &xlda);
687: if (xlda < m) {
688: MatDenseGetArrayRead(B,&b);
689: PetscMalloc1(nrhs * m, &y);
690: if (blda == m) {
691: PetscArraycpy(y,b,blda*nrhs);
692: } else {
693: for (PetscInt j = 0; j < nrhs; j++) {
694: PetscArraycpy(&y[j*m],&b[j*blda],m);
695: }
696: }
697: ylda = m;
698: MatDenseRestoreArrayRead(B,&b);
699: } else {
700: if (blda == xlda) {
701: MatCopy(B, X, SAME_NONZERO_PATTERN);
702: MatDenseGetArray(X,&y);
703: } else {
704: MatDenseGetArray(X,&y);
705: MatDenseGetArrayRead(B,&b);
706: for (PetscInt j = 0; j < nrhs; j++) {
707: PetscArraycpy(&y[j*xlda],&b[j*blda],m);
708: }
709: MatDenseRestoreArrayRead(B,&b);
710: }
711: ylda = xlda;
712: }
713: *_y = y;
714: *_ylda = ylda;
715: *_k = k;
716: *_m = m;
717: *_nrhs = nrhs;
718: return(0);
719: }
721: static PetscErrorCode MatMatSolve_SeqDense_TearDown(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ylda, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
722: {
723: PetscScalar *y;
724: PetscInt _xlda;
725: PetscBLASInt k,ylda,nrhs,xlda=0;
726: PetscErrorCode ierr;
729: y = *_y;
730: *_y = NULL;
731: k = *_k;
732: ylda = *_ylda;
733: nrhs = *_nrhs;
734: MatDenseGetLDA(X,&_xlda);
735: PetscBLASIntCast(_xlda, &xlda);
736: if (xlda != ylda) {
737: PetscScalar *xv;
738: MatDenseGetArray(X,&xv);
739: for (PetscInt j = 0; j < nrhs; j++) {
740: PetscArraycpy(&xv[j*xlda],&y[j*ylda],k);
741: }
742: MatDenseRestoreArray(X,&xv);
743: PetscFree(y);
744: } else {
745: MatDenseRestoreArray(X,&y);
746: }
747: return(0);
748: }
750: static PetscErrorCode MatMatSolve_SeqDense_LU(Mat A, Mat B, Mat X)
751: {
752: PetscScalar *y;
753: PetscBLASInt m, k, ylda, nrhs;
757: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ylda, &m, &nrhs, &k);
758: MatSolve_SeqDense_Internal_LU(A, y, ylda, m, nrhs, k, PETSC_FALSE);
759: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ylda, &m, &nrhs, &k);
760: return(0);
761: }
763: static PetscErrorCode MatMatSolveTranspose_SeqDense_LU(Mat A, Mat B, Mat X)
764: {
765: PetscScalar *y;
766: PetscBLASInt m, k, ylda, nrhs;
770: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ylda, &m, &nrhs, &k);
771: MatSolve_SeqDense_Internal_LU(A, y, ylda, m, nrhs, k, PETSC_TRUE);
772: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ylda, &m, &nrhs, &k);
773: return(0);
774: }
776: static PetscErrorCode MatMatSolve_SeqDense_Cholesky(Mat A, Mat B, Mat X)
777: {
778: PetscScalar *y;
779: PetscBLASInt m, k, ylda, nrhs;
783: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ylda, &m, &nrhs, &k);
784: MatSolve_SeqDense_Internal_Cholesky(A, y, ylda, m, nrhs, k, PETSC_FALSE);
785: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ylda, &m, &nrhs, &k);
786: return(0);
787: }
789: static PetscErrorCode MatMatSolveTranspose_SeqDense_Cholesky(Mat A, Mat B, Mat X)
790: {
791: PetscScalar *y;
792: PetscBLASInt m, k, ylda, nrhs;
796: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ylda, &m, &nrhs, &k);
797: MatSolve_SeqDense_Internal_Cholesky(A, y, ylda, m, nrhs, k, PETSC_TRUE);
798: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ylda, &m, &nrhs, &k);
799: return(0);
800: }
802: static PetscErrorCode MatMatSolve_SeqDense_QR(Mat A, Mat B, Mat X)
803: {
804: PetscScalar *y;
805: PetscBLASInt m, k, ylda, nrhs;
809: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ylda, &m, &nrhs, &k);
810: MatSolve_SeqDense_Internal_QR(A, y, ylda, m, nrhs, k);
811: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ylda, &m, &nrhs, &k);
812: return(0);
813: }
815: static PetscErrorCode MatMatSolveTranspose_SeqDense_QR(Mat A, Mat B, Mat X)
816: {
817: PetscScalar *y;
818: PetscBLASInt m, k, ylda, nrhs;
822: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ylda, &m, &nrhs, &k);
823: MatSolveTranspose_SeqDense_Internal_QR(A, y, ylda, m, nrhs, k);
824: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ylda, &m, &nrhs, &k);
825: return(0);
826: }
828: static PetscErrorCode MatConjugate_SeqDense(Mat);
830: /* ---------------------------------------------------------------*/
831: /* COMMENT: I have chosen to hide row permutation in the pivots,
832: rather than put it in the Mat->row slot.*/
833: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
834: {
835: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
837: PetscBLASInt n,m,info;
840: PetscBLASIntCast(A->cmap->n,&n);
841: PetscBLASIntCast(A->rmap->n,&m);
842: if (!mat->pivots) {
843: PetscMalloc1(A->rmap->n,&mat->pivots);
844: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
845: }
846: if (!A->rmap->n || !A->cmap->n) return(0);
847: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
848: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
849: PetscFPTrapPop();
851: if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
852: if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
854: A->ops->solve = MatSolve_SeqDense_LU;
855: A->ops->matsolve = MatMatSolve_SeqDense_LU;
856: A->ops->solvetranspose = MatSolveTranspose_SeqDense_LU;
857: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_LU;
858: A->factortype = MAT_FACTOR_LU;
860: PetscFree(A->solvertype);
861: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
863: PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
864: return(0);
865: }
867: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
868: {
869: MatFactorInfo info;
873: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
874: (*fact->ops->lufactor)(fact,NULL,NULL,&info);
875: return(0);
876: }
878: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
879: {
881: fact->preallocated = PETSC_TRUE;
882: fact->assembled = PETSC_TRUE;
883: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
884: fact->ops->solve = MatSolve_SeqDense_LU;
885: fact->ops->matsolve = MatMatSolve_SeqDense_LU;
886: fact->ops->solvetranspose = MatSolveTranspose_SeqDense_LU;
887: return(0);
888: }
890: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
891: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
892: {
893: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
895: PetscBLASInt info,n;
898: PetscBLASIntCast(A->cmap->n,&n);
899: if (!A->rmap->n || !A->cmap->n) return(0);
900: if (A->spd) {
901: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
902: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
903: PetscFPTrapPop();
904: #if defined(PETSC_USE_COMPLEX)
905: } else if (A->hermitian) {
906: if (!mat->pivots) {
907: PetscMalloc1(A->rmap->n,&mat->pivots);
908: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
909: }
910: if (!mat->fwork) {
911: PetscScalar dummy;
913: mat->lfwork = -1;
914: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
915: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
916: PetscFPTrapPop();
917: mat->lfwork = (PetscInt)PetscRealPart(dummy);
918: PetscMalloc1(mat->lfwork,&mat->fwork);
919: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
920: }
921: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
922: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
923: PetscFPTrapPop();
924: #endif
925: } else { /* symmetric case */
926: if (!mat->pivots) {
927: PetscMalloc1(A->rmap->n,&mat->pivots);
928: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
929: }
930: if (!mat->fwork) {
931: PetscScalar dummy;
933: mat->lfwork = -1;
934: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
935: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
936: PetscFPTrapPop();
937: mat->lfwork = (PetscInt)PetscRealPart(dummy);
938: PetscMalloc1(mat->lfwork,&mat->fwork);
939: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
940: }
941: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
942: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
943: PetscFPTrapPop();
944: }
945: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
947: A->ops->solve = MatSolve_SeqDense_Cholesky;
948: A->ops->matsolve = MatMatSolve_SeqDense_Cholesky;
949: A->ops->solvetranspose = MatSolveTranspose_SeqDense_Cholesky;
950: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky;
951: A->factortype = MAT_FACTOR_CHOLESKY;
953: PetscFree(A->solvertype);
954: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
956: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
957: return(0);
958: }
960: static PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
961: {
963: MatFactorInfo info;
966: info.fill = 1.0;
968: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
969: (*fact->ops->choleskyfactor)(fact,NULL,&info);
970: return(0);
971: }
973: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
974: {
976: fact->assembled = PETSC_TRUE;
977: fact->preallocated = PETSC_TRUE;
978: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
979: fact->ops->solve = MatSolve_SeqDense_Cholesky;
980: fact->ops->matsolve = MatMatSolve_SeqDense_Cholesky;
981: fact->ops->solvetranspose = MatSolve_SeqDense_Cholesky;
982: return(0);
983: }
985: static PetscErrorCode MatQRFactor_SeqDense(Mat A,IS col,const MatFactorInfo *minfo)
986: {
987: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
989: PetscBLASInt n,m,info, min, max;
992: PetscBLASIntCast(A->cmap->n,&n);
993: PetscBLASIntCast(A->rmap->n,&m);
994: max = PetscMax(m, n);
995: min = PetscMin(m, n);
996: if (!mat->tau) {
997: PetscMalloc1(min,&mat->tau);
998: PetscLogObjectMemory((PetscObject)A,min*sizeof(PetscScalar));
999: }
1000: if (!mat->pivots) {
1001: PetscMalloc1(m,&mat->pivots);
1002: PetscLogObjectMemory((PetscObject)A,m*sizeof(PetscScalar));
1003: }
1004: if (!mat->qrrhs) {
1005: MatCreateVecs(A, NULL, &(mat->qrrhs));
1006: }
1007: if (!A->rmap->n || !A->cmap->n) return(0);
1008: if (!mat->fwork) {
1009: PetscScalar dummy;
1011: mat->lfwork = -1;
1012: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1013: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,&dummy,&mat->lfwork,&info));
1014: PetscFPTrapPop();
1015: mat->lfwork = (PetscInt)PetscRealPart(dummy);
1016: PetscMalloc1(mat->lfwork,&mat->fwork);
1017: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
1018: }
1019: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1020: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,mat->fwork,&mat->lfwork,&info));
1021: PetscFPTrapPop();
1022: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to QR factorization");
1023: // TODO: try to estimate rank or test for and use geqp3 for rank revealing QR. For now just say rank is min of m and n
1024: mat->rank = min;
1026: A->ops->solve = MatSolve_SeqDense_QR;
1027: A->ops->matsolve = MatMatSolve_SeqDense_QR;
1028: A->factortype = MAT_FACTOR_QR;
1029: if (m == n) {
1030: A->ops->solvetranspose = MatSolveTranspose_SeqDense_QR;
1031: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
1032: }
1034: PetscFree(A->solvertype);
1035: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
1037: PetscLogFlops(2.0*min*min*(max-min/3.0));
1038: return(0);
1039: }
1041: static PetscErrorCode MatQRFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
1042: {
1044: MatFactorInfo info;
1047: info.fill = 1.0;
1049: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
1050: MatQRFactor_SeqDense(fact,NULL,&info);
1051: return(0);
1052: }
1054: static PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
1055: {
1059: fact->assembled = PETSC_TRUE;
1060: fact->preallocated = PETSC_TRUE;
1061: fact->ops->solve = MatSolve_SeqDense_QR;
1062: fact->ops->matsolve = MatMatSolve_SeqDense_QR;
1063: if (A->cmap->n == A->rmap->n) {
1064: fact->ops->solvetranspose = MatSolveTranspose_SeqDense_QR;
1065: fact->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
1066: }
1067: PetscObjectComposeFunction((PetscObject)fact,"MatQRFactorNumeric_C",MatQRFactorNumeric_SeqDense);
1068: return(0);
1069: }
1072: /* uses LAPACK */
1073: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
1074: {
1078: MatCreate(PetscObjectComm((PetscObject)A),fact);
1079: MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
1080: MatSetType(*fact,MATDENSE);
1081: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
1082: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
1083: (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
1084: } else {
1085: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
1086: }
1087: (*fact)->factortype = ftype;
1089: PetscFree((*fact)->solvertype);
1090: PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
1091: return(0);
1092: }
1094: /* ------------------------------------------------------------------*/
1095: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
1096: {
1097: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1098: PetscScalar *x,*v = mat->v,zero = 0.0,xt;
1099: const PetscScalar *b;
1100: PetscErrorCode ierr;
1101: PetscInt m = A->rmap->n,i;
1102: PetscBLASInt o = 1,bm = 0;
1105: #if defined(PETSC_HAVE_CUDA)
1106: if (A->offloadmask == PETSC_OFFLOAD_GPU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented");
1107: #endif
1108: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
1109: PetscBLASIntCast(m,&bm);
1110: if (flag & SOR_ZERO_INITIAL_GUESS) {
1111: /* this is a hack fix, should have another version without the second BLASdotu */
1112: VecSet(xx,zero);
1113: }
1114: VecGetArray(xx,&x);
1115: VecGetArrayRead(bb,&b);
1116: its = its*lits;
1117: if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
1118: while (its--) {
1119: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1120: for (i=0; i<m; i++) {
1121: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
1122: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
1123: }
1124: }
1125: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1126: for (i=m-1; i>=0; i--) {
1127: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
1128: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
1129: }
1130: }
1131: }
1132: VecRestoreArrayRead(bb,&b);
1133: VecRestoreArray(xx,&x);
1134: return(0);
1135: }
1137: /* -----------------------------------------------------------------*/
1138: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
1139: {
1140: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1141: const PetscScalar *v = mat->v,*x;
1142: PetscScalar *y;
1143: PetscErrorCode ierr;
1144: PetscBLASInt m, n,_One=1;
1145: PetscScalar _DOne=1.0,_DZero=0.0;
1148: PetscBLASIntCast(A->rmap->n,&m);
1149: PetscBLASIntCast(A->cmap->n,&n);
1150: VecGetArrayRead(xx,&x);
1151: VecGetArrayWrite(yy,&y);
1152: if (!A->rmap->n || !A->cmap->n) {
1153: PetscBLASInt i;
1154: for (i=0; i<n; i++) y[i] = 0.0;
1155: } else {
1156: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
1157: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
1158: }
1159: VecRestoreArrayRead(xx,&x);
1160: VecRestoreArrayWrite(yy,&y);
1161: return(0);
1162: }
1164: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
1165: {
1166: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1167: PetscScalar *y,_DOne=1.0,_DZero=0.0;
1168: PetscErrorCode ierr;
1169: PetscBLASInt m, n, _One=1;
1170: const PetscScalar *v = mat->v,*x;
1173: PetscBLASIntCast(A->rmap->n,&m);
1174: PetscBLASIntCast(A->cmap->n,&n);
1175: VecGetArrayRead(xx,&x);
1176: VecGetArrayWrite(yy,&y);
1177: if (!A->rmap->n || !A->cmap->n) {
1178: PetscBLASInt i;
1179: for (i=0; i<m; i++) y[i] = 0.0;
1180: } else {
1181: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
1182: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
1183: }
1184: VecRestoreArrayRead(xx,&x);
1185: VecRestoreArrayWrite(yy,&y);
1186: return(0);
1187: }
1189: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
1190: {
1191: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1192: const PetscScalar *v = mat->v,*x;
1193: PetscScalar *y,_DOne=1.0;
1194: PetscErrorCode ierr;
1195: PetscBLASInt m, n, _One=1;
1198: PetscBLASIntCast(A->rmap->n,&m);
1199: PetscBLASIntCast(A->cmap->n,&n);
1200: VecCopy(zz,yy);
1201: if (!A->rmap->n || !A->cmap->n) return(0);
1202: VecGetArrayRead(xx,&x);
1203: VecGetArray(yy,&y);
1204: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
1205: VecRestoreArrayRead(xx,&x);
1206: VecRestoreArray(yy,&y);
1207: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
1208: return(0);
1209: }
1211: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
1212: {
1213: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1214: const PetscScalar *v = mat->v,*x;
1215: PetscScalar *y;
1216: PetscErrorCode ierr;
1217: PetscBLASInt m, n, _One=1;
1218: PetscScalar _DOne=1.0;
1221: PetscBLASIntCast(A->rmap->n,&m);
1222: PetscBLASIntCast(A->cmap->n,&n);
1223: VecCopy(zz,yy);
1224: if (!A->rmap->n || !A->cmap->n) return(0);
1225: VecGetArrayRead(xx,&x);
1226: VecGetArray(yy,&y);
1227: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
1228: VecRestoreArrayRead(xx,&x);
1229: VecRestoreArray(yy,&y);
1230: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
1231: return(0);
1232: }
1234: /* -----------------------------------------------------------------*/
1235: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
1236: {
1237: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1239: PetscInt i;
1242: *ncols = A->cmap->n;
1243: if (cols) {
1244: PetscMalloc1(A->cmap->n,cols);
1245: for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
1246: }
1247: if (vals) {
1248: const PetscScalar *v;
1250: MatDenseGetArrayRead(A,&v);
1251: PetscMalloc1(A->cmap->n,vals);
1252: v += row;
1253: for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
1254: MatDenseRestoreArrayRead(A,&v);
1255: }
1256: return(0);
1257: }
1259: static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
1260: {
1264: if (ncols) *ncols = 0;
1265: if (cols) {PetscFree(*cols);}
1266: if (vals) {PetscFree(*vals); }
1267: return(0);
1268: }
1269: /* ----------------------------------------------------------------*/
1270: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
1271: {
1272: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1273: PetscScalar *av;
1274: PetscInt i,j,idx=0;
1275: #if defined(PETSC_HAVE_CUDA)
1276: PetscOffloadMask oldf;
1277: #endif
1278: PetscErrorCode ierr;
1281: MatDenseGetArray(A,&av);
1282: if (!mat->roworiented) {
1283: if (addv == INSERT_VALUES) {
1284: for (j=0; j<n; j++) {
1285: if (indexn[j] < 0) {idx += m; continue;}
1286: if (PetscUnlikelyDebug(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);
1287: for (i=0; i<m; i++) {
1288: if (indexm[i] < 0) {idx++; continue;}
1289: if (PetscUnlikelyDebug(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);
1290: av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1291: }
1292: }
1293: } else {
1294: for (j=0; j<n; j++) {
1295: if (indexn[j] < 0) {idx += m; continue;}
1296: if (PetscUnlikelyDebug(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);
1297: for (i=0; i<m; i++) {
1298: if (indexm[i] < 0) {idx++; continue;}
1299: if (PetscUnlikelyDebug(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);
1300: av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1301: }
1302: }
1303: }
1304: } else {
1305: if (addv == INSERT_VALUES) {
1306: for (i=0; i<m; i++) {
1307: if (indexm[i] < 0) { idx += n; continue;}
1308: if (PetscUnlikelyDebug(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);
1309: for (j=0; j<n; j++) {
1310: if (indexn[j] < 0) { idx++; continue;}
1311: if (PetscUnlikelyDebug(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);
1312: av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1313: }
1314: }
1315: } else {
1316: for (i=0; i<m; i++) {
1317: if (indexm[i] < 0) { idx += n; continue;}
1318: if (PetscUnlikelyDebug(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);
1319: for (j=0; j<n; j++) {
1320: if (indexn[j] < 0) { idx++; continue;}
1321: if (PetscUnlikelyDebug(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);
1322: av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1323: }
1324: }
1325: }
1326: }
1327: /* hack to prevent unneeded copy to the GPU while returning the array */
1328: #if defined(PETSC_HAVE_CUDA)
1329: oldf = A->offloadmask;
1330: A->offloadmask = PETSC_OFFLOAD_GPU;
1331: #endif
1332: MatDenseRestoreArray(A,&av);
1333: #if defined(PETSC_HAVE_CUDA)
1334: A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1335: #endif
1336: return(0);
1337: }
1339: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1340: {
1341: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1342: const PetscScalar *vv;
1343: PetscInt i,j;
1344: PetscErrorCode ierr;
1347: MatDenseGetArrayRead(A,&vv);
1348: /* row-oriented output */
1349: for (i=0; i<m; i++) {
1350: if (indexm[i] < 0) {v += n;continue;}
1351: 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);
1352: for (j=0; j<n; j++) {
1353: if (indexn[j] < 0) {v++; continue;}
1354: 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);
1355: *v++ = vv[indexn[j]*mat->lda + indexm[i]];
1356: }
1357: }
1358: MatDenseRestoreArrayRead(A,&vv);
1359: return(0);
1360: }
1362: /* -----------------------------------------------------------------*/
1364: PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer)
1365: {
1366: PetscErrorCode ierr;
1367: PetscBool skipHeader;
1368: PetscViewerFormat format;
1369: PetscInt header[4],M,N,m,lda,i,j,k;
1370: const PetscScalar *v;
1371: PetscScalar *vwork;
1374: PetscViewerSetUp(viewer);
1375: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
1376: PetscViewerGetFormat(viewer,&format);
1377: if (skipHeader) format = PETSC_VIEWER_NATIVE;
1379: MatGetSize(mat,&M,&N);
1381: /* write matrix header */
1382: header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N;
1383: header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N;
1384: if (!skipHeader) {PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);}
1386: MatGetLocalSize(mat,&m,NULL);
1387: if (format != PETSC_VIEWER_NATIVE) {
1388: PetscInt nnz = m*N, *iwork;
1389: /* store row lengths for each row */
1390: PetscMalloc1(nnz,&iwork);
1391: for (i=0; i<m; i++) iwork[i] = N;
1392: PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1393: /* store column indices (zero start index) */
1394: for (k=0, i=0; i<m; i++)
1395: for (j=0; j<N; j++, k++)
1396: iwork[k] = j;
1397: PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1398: PetscFree(iwork);
1399: }
1400: /* store matrix values as a dense matrix in row major order */
1401: PetscMalloc1(m*N,&vwork);
1402: MatDenseGetArrayRead(mat,&v);
1403: MatDenseGetLDA(mat,&lda);
1404: for (k=0, i=0; i<m; i++)
1405: for (j=0; j<N; j++, k++)
1406: vwork[k] = v[i+lda*j];
1407: MatDenseRestoreArrayRead(mat,&v);
1408: PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1409: PetscFree(vwork);
1410: return(0);
1411: }
1413: PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer)
1414: {
1416: PetscBool skipHeader;
1417: PetscInt header[4],M,N,m,nz,lda,i,j,k;
1418: PetscInt rows,cols;
1419: PetscScalar *v,*vwork;
1422: PetscViewerSetUp(viewer);
1423: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
1425: if (!skipHeader) {
1426: PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
1427: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file");
1428: M = header[1]; N = header[2];
1429: if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M);
1430: if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N);
1431: nz = header[3];
1432: if (nz != MATRIX_BINARY_FORMAT_DENSE && nz < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Unknown matrix format %D in file",nz);
1433: } else {
1434: MatGetSize(mat,&M,&N);
1435: if (M < 0 || N < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Matrix binary file header was skipped, thus the user must specify the global sizes of input matrix");
1436: nz = MATRIX_BINARY_FORMAT_DENSE;
1437: }
1439: /* setup global sizes if not set */
1440: if (mat->rmap->N < 0) mat->rmap->N = M;
1441: if (mat->cmap->N < 0) mat->cmap->N = N;
1442: MatSetUp(mat);
1443: /* check if global sizes are correct */
1444: MatGetSize(mat,&rows,&cols);
1445: if (M != rows || N != cols) SETERRQ4(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols);
1447: MatGetSize(mat,NULL,&N);
1448: MatGetLocalSize(mat,&m,NULL);
1449: MatDenseGetArray(mat,&v);
1450: MatDenseGetLDA(mat,&lda);
1451: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */
1452: PetscInt nnz = m*N;
1453: /* read in matrix values */
1454: PetscMalloc1(nnz,&vwork);
1455: PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1456: /* store values in column major order */
1457: for (j=0; j<N; j++)
1458: for (i=0; i<m; i++)
1459: v[i+lda*j] = vwork[i*N+j];
1460: PetscFree(vwork);
1461: } else { /* matrix in file is sparse format */
1462: PetscInt nnz = 0, *rlens, *icols;
1463: /* read in row lengths */
1464: PetscMalloc1(m,&rlens);
1465: PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1466: for (i=0; i<m; i++) nnz += rlens[i];
1467: /* read in column indices and values */
1468: PetscMalloc2(nnz,&icols,nnz,&vwork);
1469: PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1470: PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1471: /* store values in column major order */
1472: for (k=0, i=0; i<m; i++)
1473: for (j=0; j<rlens[i]; j++, k++)
1474: v[i+lda*icols[k]] = vwork[k];
1475: PetscFree(rlens);
1476: PetscFree2(icols,vwork);
1477: }
1478: MatDenseRestoreArray(mat,&v);
1479: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
1480: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
1481: return(0);
1482: }
1484: PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1485: {
1486: PetscBool isbinary, ishdf5;
1492: /* force binary viewer to load .info file if it has not yet done so */
1493: PetscViewerSetUp(viewer);
1494: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1495: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);
1496: if (isbinary) {
1497: MatLoad_Dense_Binary(newMat,viewer);
1498: } else if (ishdf5) {
1499: #if defined(PETSC_HAVE_HDF5)
1500: MatLoad_Dense_HDF5(newMat,viewer);
1501: #else
1502: SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1503: #endif
1504: } else {
1505: 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);
1506: }
1507: return(0);
1508: }
1510: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1511: {
1512: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1513: PetscErrorCode ierr;
1514: PetscInt i,j;
1515: const char *name;
1516: PetscScalar *v,*av;
1517: PetscViewerFormat format;
1518: #if defined(PETSC_USE_COMPLEX)
1519: PetscBool allreal = PETSC_TRUE;
1520: #endif
1523: MatDenseGetArrayRead(A,(const PetscScalar**)&av);
1524: PetscViewerGetFormat(viewer,&format);
1525: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1526: return(0); /* do nothing for now */
1527: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1528: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1529: for (i=0; i<A->rmap->n; i++) {
1530: v = av + i;
1531: PetscViewerASCIIPrintf(viewer,"row %D:",i);
1532: for (j=0; j<A->cmap->n; j++) {
1533: #if defined(PETSC_USE_COMPLEX)
1534: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1535: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1536: } else if (PetscRealPart(*v)) {
1537: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1538: }
1539: #else
1540: if (*v) {
1541: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1542: }
1543: #endif
1544: v += a->lda;
1545: }
1546: PetscViewerASCIIPrintf(viewer,"\n");
1547: }
1548: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1549: } else {
1550: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1551: #if defined(PETSC_USE_COMPLEX)
1552: /* determine if matrix has all real values */
1553: v = av;
1554: for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1555: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1556: }
1557: #endif
1558: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1559: PetscObjectGetName((PetscObject)A,&name);
1560: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1561: PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1562: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1563: }
1565: for (i=0; i<A->rmap->n; i++) {
1566: v = av + i;
1567: for (j=0; j<A->cmap->n; j++) {
1568: #if defined(PETSC_USE_COMPLEX)
1569: if (allreal) {
1570: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1571: } else {
1572: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1573: }
1574: #else
1575: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1576: #endif
1577: v += a->lda;
1578: }
1579: PetscViewerASCIIPrintf(viewer,"\n");
1580: }
1581: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1582: PetscViewerASCIIPrintf(viewer,"];\n");
1583: }
1584: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1585: }
1586: MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);
1587: PetscViewerFlush(viewer);
1588: return(0);
1589: }
1591: #include <petscdraw.h>
1592: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1593: {
1594: Mat A = (Mat) Aa;
1595: PetscErrorCode ierr;
1596: PetscInt m = A->rmap->n,n = A->cmap->n,i,j;
1597: int color = PETSC_DRAW_WHITE;
1598: const PetscScalar *v;
1599: PetscViewer viewer;
1600: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1601: PetscViewerFormat format;
1604: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1605: PetscViewerGetFormat(viewer,&format);
1606: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
1608: /* Loop over matrix elements drawing boxes */
1609: MatDenseGetArrayRead(A,&v);
1610: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1611: PetscDrawCollectiveBegin(draw);
1612: /* Blue for negative and Red for positive */
1613: for (j = 0; j < n; j++) {
1614: x_l = j; x_r = x_l + 1.0;
1615: for (i = 0; i < m; i++) {
1616: y_l = m - i - 1.0;
1617: y_r = y_l + 1.0;
1618: if (PetscRealPart(v[j*m+i]) > 0.) color = PETSC_DRAW_RED;
1619: else if (PetscRealPart(v[j*m+i]) < 0.) color = PETSC_DRAW_BLUE;
1620: else continue;
1621: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1622: }
1623: }
1624: PetscDrawCollectiveEnd(draw);
1625: } else {
1626: /* use contour shading to indicate magnitude of values */
1627: /* first determine max of all nonzero values */
1628: PetscReal minv = 0.0, maxv = 0.0;
1629: PetscDraw popup;
1631: for (i=0; i < m*n; i++) {
1632: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1633: }
1634: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1635: PetscDrawGetPopup(draw,&popup);
1636: PetscDrawScalePopup(popup,minv,maxv);
1638: PetscDrawCollectiveBegin(draw);
1639: for (j=0; j<n; j++) {
1640: x_l = j;
1641: x_r = x_l + 1.0;
1642: for (i=0; i<m; i++) {
1643: y_l = m - i - 1.0;
1644: y_r = y_l + 1.0;
1645: color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1646: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1647: }
1648: }
1649: PetscDrawCollectiveEnd(draw);
1650: }
1651: MatDenseRestoreArrayRead(A,&v);
1652: return(0);
1653: }
1655: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1656: {
1657: PetscDraw draw;
1658: PetscBool isnull;
1659: PetscReal xr,yr,xl,yl,h,w;
1663: PetscViewerDrawGetDraw(viewer,0,&draw);
1664: PetscDrawIsNull(draw,&isnull);
1665: if (isnull) return(0);
1667: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1668: xr += w; yr += h; xl = -w; yl = -h;
1669: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1670: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1671: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1672: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1673: PetscDrawSave(draw);
1674: return(0);
1675: }
1677: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1678: {
1680: PetscBool iascii,isbinary,isdraw;
1683: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1684: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1685: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1687: if (iascii) {
1688: MatView_SeqDense_ASCII(A,viewer);
1689: } else if (isbinary) {
1690: MatView_Dense_Binary(A,viewer);
1691: } else if (isdraw) {
1692: MatView_SeqDense_Draw(A,viewer);
1693: }
1694: return(0);
1695: }
1697: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar *array)
1698: {
1699: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1702: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1703: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1704: if (a->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreArray() first");
1705: a->unplacedarray = a->v;
1706: a->unplaced_user_alloc = a->user_alloc;
1707: a->v = (PetscScalar*) array;
1708: a->user_alloc = PETSC_TRUE;
1709: #if defined(PETSC_HAVE_CUDA)
1710: A->offloadmask = PETSC_OFFLOAD_CPU;
1711: #endif
1712: return(0);
1713: }
1715: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1716: {
1717: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1720: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1721: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1722: a->v = a->unplacedarray;
1723: a->user_alloc = a->unplaced_user_alloc;
1724: a->unplacedarray = NULL;
1725: #if defined(PETSC_HAVE_CUDA)
1726: A->offloadmask = PETSC_OFFLOAD_CPU;
1727: #endif
1728: return(0);
1729: }
1731: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A,const PetscScalar *array)
1732: {
1733: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1737: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1738: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1739: if (!a->user_alloc) { PetscFree(a->v); }
1740: a->v = (PetscScalar*) array;
1741: a->user_alloc = PETSC_FALSE;
1742: #if defined(PETSC_HAVE_CUDA)
1743: A->offloadmask = PETSC_OFFLOAD_CPU;
1744: #endif
1745: return(0);
1746: }
1748: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1749: {
1750: Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
1754: #if defined(PETSC_USE_LOG)
1755: PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1756: #endif
1757: VecDestroy(&(l->qrrhs));
1758: PetscFree(l->tau);
1759: PetscFree(l->pivots);
1760: PetscFree(l->fwork);
1761: MatDestroy(&l->ptapwork);
1762: if (!l->user_alloc) {PetscFree(l->v);}
1763: if (!l->unplaced_user_alloc) {PetscFree(l->unplacedarray);}
1764: if (l->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1765: if (l->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1766: VecDestroy(&l->cvec);
1767: MatDestroy(&l->cmat);
1768: PetscFree(mat->data);
1770: PetscObjectChangeTypeName((PetscObject)mat,NULL);
1771: PetscObjectComposeFunction((PetscObject)mat,"MatQRFactor_C",NULL);
1772: PetscObjectComposeFunction((PetscObject)mat,"MatQRFactorNumeric_C",NULL);
1773: PetscObjectComposeFunction((PetscObject)mat,"MatQRFactorSymbolic_C",NULL);
1774: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);
1775: PetscObjectComposeFunction((PetscObject)mat,"MatDenseSetLDA_C",NULL);
1776: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1777: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1778: PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1779: PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1780: PetscObjectComposeFunction((PetscObject)mat,"MatDenseReplaceArray_C",NULL);
1781: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);
1782: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);
1783: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayWrite_C",NULL);
1784: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayWrite_C",NULL);
1785: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1786: #if defined(PETSC_HAVE_ELEMENTAL)
1787: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1788: #endif
1789: #if defined(PETSC_HAVE_SCALAPACK)
1790: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_scalapack_C",NULL);
1791: #endif
1792: #if defined(PETSC_HAVE_CUDA)
1793: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);
1794: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL);
1795: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL);
1796: #endif
1797: PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1798: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL);
1799: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL);
1800: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL);
1801: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL);
1803: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
1804: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
1805: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVec_C",NULL);
1806: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVec_C",NULL);
1807: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecRead_C",NULL);
1808: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecRead_C",NULL);
1809: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecWrite_C",NULL);
1810: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecWrite_C",NULL);
1811: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetSubMatrix_C",NULL);
1812: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreSubMatrix_C",NULL);
1813: return(0);
1814: }
1816: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1817: {
1818: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1820: PetscInt k,j,m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1821: PetscScalar *v,tmp;
1824: if (reuse == MAT_INPLACE_MATRIX) {
1825: if (m == n) { /* in place transpose */
1826: MatDenseGetArray(A,&v);
1827: for (j=0; j<m; j++) {
1828: for (k=0; k<j; k++) {
1829: tmp = v[j + k*M];
1830: v[j + k*M] = v[k + j*M];
1831: v[k + j*M] = tmp;
1832: }
1833: }
1834: MatDenseRestoreArray(A,&v);
1835: } else { /* reuse memory, temporary allocates new memory */
1836: PetscScalar *v2;
1837: PetscLayout tmplayout;
1839: PetscMalloc1((size_t)m*n,&v2);
1840: MatDenseGetArray(A,&v);
1841: for (j=0; j<n; j++) {
1842: for (k=0; k<m; k++) v2[j + (size_t)k*n] = v[k + (size_t)j*M];
1843: }
1844: PetscArraycpy(v,v2,(size_t)m*n);
1845: PetscFree(v2);
1846: MatDenseRestoreArray(A,&v);
1847: /* cleanup size dependent quantities */
1848: VecDestroy(&mat->cvec);
1849: MatDestroy(&mat->cmat);
1850: PetscFree(mat->pivots);
1851: PetscFree(mat->fwork);
1852: MatDestroy(&mat->ptapwork);
1853: /* swap row/col layouts */
1854: mat->lda = n;
1855: tmplayout = A->rmap;
1856: A->rmap = A->cmap;
1857: A->cmap = tmplayout;
1858: }
1859: } else { /* out-of-place transpose */
1860: Mat tmat;
1861: Mat_SeqDense *tmatd;
1862: PetscScalar *v2;
1863: PetscInt M2;
1865: if (reuse == MAT_INITIAL_MATRIX) {
1866: MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1867: MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1868: MatSetType(tmat,((PetscObject)A)->type_name);
1869: MatSeqDenseSetPreallocation(tmat,NULL);
1870: } else tmat = *matout;
1872: MatDenseGetArrayRead(A,(const PetscScalar**)&v);
1873: MatDenseGetArray(tmat,&v2);
1874: tmatd = (Mat_SeqDense*)tmat->data;
1875: M2 = tmatd->lda;
1876: for (j=0; j<n; j++) {
1877: for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1878: }
1879: MatDenseRestoreArray(tmat,&v2);
1880: MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);
1881: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1882: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1883: *matout = tmat;
1884: }
1885: return(0);
1886: }
1888: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg)
1889: {
1890: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1891: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1892: PetscInt i;
1893: const PetscScalar *v1,*v2;
1894: PetscErrorCode ierr;
1897: if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1898: if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1899: MatDenseGetArrayRead(A1,&v1);
1900: MatDenseGetArrayRead(A2,&v2);
1901: for (i=0; i<A1->cmap->n; i++) {
1902: PetscArraycmp(v1,v2,A1->rmap->n,flg);
1903: if (*flg == PETSC_FALSE) return(0);
1904: v1 += mat1->lda;
1905: v2 += mat2->lda;
1906: }
1907: MatDenseRestoreArrayRead(A1,&v1);
1908: MatDenseRestoreArrayRead(A2,&v2);
1909: *flg = PETSC_TRUE;
1910: return(0);
1911: }
1913: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1914: {
1915: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1916: PetscInt i,n,len;
1917: PetscScalar *x;
1918: const PetscScalar *vv;
1919: PetscErrorCode ierr;
1922: VecGetSize(v,&n);
1923: VecGetArray(v,&x);
1924: len = PetscMin(A->rmap->n,A->cmap->n);
1925: MatDenseGetArrayRead(A,&vv);
1926: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1927: for (i=0; i<len; i++) {
1928: x[i] = vv[i*mat->lda + i];
1929: }
1930: MatDenseRestoreArrayRead(A,&vv);
1931: VecRestoreArray(v,&x);
1932: return(0);
1933: }
1935: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1936: {
1937: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1938: const PetscScalar *l,*r;
1939: PetscScalar x,*v,*vv;
1940: PetscErrorCode ierr;
1941: PetscInt i,j,m = A->rmap->n,n = A->cmap->n;
1944: MatDenseGetArray(A,&vv);
1945: if (ll) {
1946: VecGetSize(ll,&m);
1947: VecGetArrayRead(ll,&l);
1948: if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1949: for (i=0; i<m; i++) {
1950: x = l[i];
1951: v = vv + i;
1952: for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1953: }
1954: VecRestoreArrayRead(ll,&l);
1955: PetscLogFlops(1.0*n*m);
1956: }
1957: if (rr) {
1958: VecGetSize(rr,&n);
1959: VecGetArrayRead(rr,&r);
1960: if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1961: for (i=0; i<n; i++) {
1962: x = r[i];
1963: v = vv + i*mat->lda;
1964: for (j=0; j<m; j++) (*v++) *= x;
1965: }
1966: VecRestoreArrayRead(rr,&r);
1967: PetscLogFlops(1.0*n*m);
1968: }
1969: MatDenseRestoreArray(A,&vv);
1970: return(0);
1971: }
1973: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1974: {
1975: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1976: PetscScalar *v,*vv;
1977: PetscReal sum = 0.0;
1978: PetscInt lda =mat->lda,m=A->rmap->n,i,j;
1979: PetscErrorCode ierr;
1982: MatDenseGetArrayRead(A,(const PetscScalar**)&vv);
1983: v = vv;
1984: if (type == NORM_FROBENIUS) {
1985: if (lda>m) {
1986: for (j=0; j<A->cmap->n; j++) {
1987: v = vv+j*lda;
1988: for (i=0; i<m; i++) {
1989: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1990: }
1991: }
1992: } else {
1993: #if defined(PETSC_USE_REAL___FP16)
1994: PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1995: PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&cnt,v,&one));
1996: }
1997: #else
1998: for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1999: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
2000: }
2001: }
2002: *nrm = PetscSqrtReal(sum);
2003: #endif
2004: PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
2005: } else if (type == NORM_1) {
2006: *nrm = 0.0;
2007: for (j=0; j<A->cmap->n; j++) {
2008: v = vv + j*mat->lda;
2009: sum = 0.0;
2010: for (i=0; i<A->rmap->n; i++) {
2011: sum += PetscAbsScalar(*v); v++;
2012: }
2013: if (sum > *nrm) *nrm = sum;
2014: }
2015: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
2016: } else if (type == NORM_INFINITY) {
2017: *nrm = 0.0;
2018: for (j=0; j<A->rmap->n; j++) {
2019: v = vv + j;
2020: sum = 0.0;
2021: for (i=0; i<A->cmap->n; i++) {
2022: sum += PetscAbsScalar(*v); v += mat->lda;
2023: }
2024: if (sum > *nrm) *nrm = sum;
2025: }
2026: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
2027: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
2028: MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);
2029: return(0);
2030: }
2032: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
2033: {
2034: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
2038: switch (op) {
2039: case MAT_ROW_ORIENTED:
2040: aij->roworiented = flg;
2041: break;
2042: case MAT_NEW_NONZERO_LOCATIONS:
2043: case MAT_NEW_NONZERO_LOCATION_ERR:
2044: case MAT_NEW_NONZERO_ALLOCATION_ERR:
2045: case MAT_FORCE_DIAGONAL_ENTRIES:
2046: case MAT_KEEP_NONZERO_PATTERN:
2047: case MAT_IGNORE_OFF_PROC_ENTRIES:
2048: case MAT_USE_HASH_TABLE:
2049: case MAT_IGNORE_ZERO_ENTRIES:
2050: case MAT_IGNORE_LOWER_TRIANGULAR:
2051: case MAT_SORTED_FULL:
2052: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
2053: break;
2054: case MAT_SPD:
2055: case MAT_SYMMETRIC:
2056: case MAT_STRUCTURALLY_SYMMETRIC:
2057: case MAT_HERMITIAN:
2058: case MAT_SYMMETRY_ETERNAL:
2059: /* These options are handled directly by MatSetOption() */
2060: break;
2061: default:
2062: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
2063: }
2064: return(0);
2065: }
2067: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
2068: {
2069: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
2071: PetscInt lda=l->lda,m=A->rmap->n,j;
2072: PetscScalar *v;
2075: MatDenseGetArray(A,&v);
2076: if (lda>m) {
2077: for (j=0; j<A->cmap->n; j++) {
2078: PetscArrayzero(v+j*lda,m);
2079: }
2080: } else {
2081: PetscArrayzero(v,A->rmap->n*A->cmap->n);
2082: }
2083: MatDenseRestoreArray(A,&v);
2084: return(0);
2085: }
2087: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2088: {
2089: PetscErrorCode ierr;
2090: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
2091: PetscInt m = l->lda, n = A->cmap->n, i,j;
2092: PetscScalar *slot,*bb,*v;
2093: const PetscScalar *xx;
2096: if (PetscDefined(USE_DEBUG)) {
2097: for (i=0; i<N; i++) {
2098: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
2099: 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);
2100: }
2101: }
2102: if (!N) return(0);
2104: /* fix right hand side if needed */
2105: if (x && b) {
2106: VecGetArrayRead(x,&xx);
2107: VecGetArray(b,&bb);
2108: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
2109: VecRestoreArrayRead(x,&xx);
2110: VecRestoreArray(b,&bb);
2111: }
2113: MatDenseGetArray(A,&v);
2114: for (i=0; i<N; i++) {
2115: slot = v + rows[i];
2116: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
2117: }
2118: if (diag != 0.0) {
2119: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
2120: for (i=0; i<N; i++) {
2121: slot = v + (m+1)*rows[i];
2122: *slot = diag;
2123: }
2124: }
2125: MatDenseRestoreArray(A,&v);
2126: return(0);
2127: }
2129: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda)
2130: {
2131: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
2134: *lda = mat->lda;
2135: return(0);
2136: }
2138: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array)
2139: {
2140: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
2143: if (mat->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
2144: *array = mat->v;
2145: return(0);
2146: }
2148: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array)
2149: {
2151: *array = NULL;
2152: return(0);
2153: }
2155: /*@C
2156: MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray()
2158: Not collective
2160: Input Parameter:
2161: . mat - a MATSEQDENSE or MATMPIDENSE matrix
2163: Output Parameter:
2164: . lda - the leading dimension
2166: Level: intermediate
2168: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseSetLDA()
2169: @*/
2170: PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda)
2171: {
2177: PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));
2178: return(0);
2179: }
2181: /*@C
2182: MatDenseSetLDA - Sets the leading dimension of the array used by the dense matrix
2184: Not collective
2186: Input Parameter:
2187: + mat - a MATSEQDENSE or MATMPIDENSE matrix
2188: - lda - the leading dimension
2190: Level: intermediate
2192: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetLDA()
2193: @*/
2194: PetscErrorCode MatDenseSetLDA(Mat A,PetscInt lda)
2195: {
2200: PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda));
2201: return(0);
2202: }
2204: /*@C
2205: MatDenseGetArray - gives read-write access to the array where the data for a dense matrix is stored
2207: Logically Collective on Mat
2209: Input Parameter:
2210: . mat - a dense matrix
2212: Output Parameter:
2213: . array - pointer to the data
2215: Level: intermediate
2217: .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2218: @*/
2219: PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array)
2220: {
2226: PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
2227: return(0);
2228: }
2230: /*@C
2231: MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
2233: Logically Collective on Mat
2235: Input Parameters:
2236: + mat - a dense matrix
2237: - array - pointer to the data
2239: Level: intermediate
2241: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2242: @*/
2243: PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array)
2244: {
2250: PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
2251: PetscObjectStateIncrease((PetscObject)A);
2252: #if defined(PETSC_HAVE_CUDA)
2253: A->offloadmask = PETSC_OFFLOAD_CPU;
2254: #endif
2255: return(0);
2256: }
2258: /*@C
2259: MatDenseGetArrayRead - gives read-only access to the array where the data for a dense matrix is stored
2261: Not Collective
2263: Input Parameter:
2264: . mat - a dense matrix
2266: Output Parameter:
2267: . array - pointer to the data
2269: Level: intermediate
2271: .seealso: MatDenseRestoreArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2272: @*/
2273: PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array)
2274: {
2280: PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
2281: return(0);
2282: }
2284: /*@C
2285: MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayRead()
2287: Not Collective
2289: Input Parameters:
2290: + mat - a dense matrix
2291: - array - pointer to the data
2293: Level: intermediate
2295: .seealso: MatDenseGetArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2296: @*/
2297: PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
2298: {
2304: PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
2305: return(0);
2306: }
2308: /*@C
2309: MatDenseGetArrayWrite - gives write-only access to the array where the data for a dense matrix is stored
2311: Not Collective
2313: Input Parameter:
2314: . mat - a dense matrix
2316: Output Parameter:
2317: . array - pointer to the data
2319: Level: intermediate
2321: .seealso: MatDenseRestoreArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
2322: @*/
2323: PetscErrorCode MatDenseGetArrayWrite(Mat A,PetscScalar **array)
2324: {
2330: PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array));
2331: return(0);
2332: }
2334: /*@C
2335: MatDenseRestoreArrayWrite - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayWrite()
2337: Not Collective
2339: Input Parameters:
2340: + mat - a dense matrix
2341: - array - pointer to the data
2343: Level: intermediate
2345: .seealso: MatDenseGetArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
2346: @*/
2347: PetscErrorCode MatDenseRestoreArrayWrite(Mat A,PetscScalar **array)
2348: {
2354: PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array));
2355: PetscObjectStateIncrease((PetscObject)A);
2356: #if defined(PETSC_HAVE_CUDA)
2357: A->offloadmask = PETSC_OFFLOAD_CPU;
2358: #endif
2359: return(0);
2360: }
2362: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B)
2363: {
2364: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
2366: PetscInt i,j,nrows,ncols,blda;
2367: const PetscInt *irow,*icol;
2368: PetscScalar *av,*bv,*v = mat->v;
2369: Mat newmat;
2372: ISGetIndices(isrow,&irow);
2373: ISGetIndices(iscol,&icol);
2374: ISGetLocalSize(isrow,&nrows);
2375: ISGetLocalSize(iscol,&ncols);
2377: /* Check submatrixcall */
2378: if (scall == MAT_REUSE_MATRIX) {
2379: PetscInt n_cols,n_rows;
2380: MatGetSize(*B,&n_rows,&n_cols);
2381: if (n_rows != nrows || n_cols != ncols) {
2382: /* resize the result matrix to match number of requested rows/columns */
2383: MatSetSizes(*B,nrows,ncols,nrows,ncols);
2384: }
2385: newmat = *B;
2386: } else {
2387: /* Create and fill new matrix */
2388: MatCreate(PetscObjectComm((PetscObject)A),&newmat);
2389: MatSetSizes(newmat,nrows,ncols,nrows,ncols);
2390: MatSetType(newmat,((PetscObject)A)->type_name);
2391: MatSeqDenseSetPreallocation(newmat,NULL);
2392: }
2394: /* Now extract the data pointers and do the copy,column at a time */
2395: MatDenseGetArray(newmat,&bv);
2396: MatDenseGetLDA(newmat,&blda);
2397: for (i=0; i<ncols; i++) {
2398: av = v + mat->lda*icol[i];
2399: for (j=0; j<nrows; j++) bv[j] = av[irow[j]];
2400: bv += blda;
2401: }
2402: MatDenseRestoreArray(newmat,&bv);
2404: /* Assemble the matrices so that the correct flags are set */
2405: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
2406: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
2408: /* Free work space */
2409: ISRestoreIndices(isrow,&irow);
2410: ISRestoreIndices(iscol,&icol);
2411: *B = newmat;
2412: return(0);
2413: }
2415: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
2416: {
2418: PetscInt i;
2421: if (scall == MAT_INITIAL_MATRIX) {
2422: PetscCalloc1(n,B);
2423: }
2425: for (i=0; i<n; i++) {
2426: MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],scall,&(*B)[i]);
2427: }
2428: return(0);
2429: }
2431: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
2432: {
2434: return(0);
2435: }
2437: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
2438: {
2440: return(0);
2441: }
2443: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
2444: {
2445: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
2446: PetscErrorCode ierr;
2447: const PetscScalar *va;
2448: PetscScalar *vb;
2449: PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;
2452: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2453: if (A->ops->copy != B->ops->copy) {
2454: MatCopy_Basic(A,B,str);
2455: return(0);
2456: }
2457: if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
2458: MatDenseGetArrayRead(A,&va);
2459: MatDenseGetArray(B,&vb);
2460: if (lda1>m || lda2>m) {
2461: for (j=0; j<n; j++) {
2462: PetscArraycpy(vb+j*lda2,va+j*lda1,m);
2463: }
2464: } else {
2465: PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);
2466: }
2467: MatDenseRestoreArray(B,&vb);
2468: MatDenseRestoreArrayRead(A,&va);
2469: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2470: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2471: return(0);
2472: }
2474: static PetscErrorCode MatSetUp_SeqDense(Mat A)
2475: {
2479: PetscLayoutSetUp(A->rmap);
2480: PetscLayoutSetUp(A->cmap);
2481: if (!A->preallocated) {
2482: MatSeqDenseSetPreallocation(A,NULL);
2483: }
2484: return(0);
2485: }
2487: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2488: {
2489: Mat_SeqDense *mat = (Mat_SeqDense *) A->data;
2490: PetscInt i,nz = A->rmap->n*A->cmap->n;
2491: PetscInt min = PetscMin(A->rmap->n,A->cmap->n);
2492: PetscScalar *aa;
2496: MatDenseGetArray(A,&aa);
2497: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
2498: MatDenseRestoreArray(A,&aa);
2499: if (mat->tau) for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]);
2500: return(0);
2501: }
2503: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2504: {
2505: PetscInt i,nz = A->rmap->n*A->cmap->n;
2506: PetscScalar *aa;
2510: MatDenseGetArray(A,&aa);
2511: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
2512: MatDenseRestoreArray(A,&aa);
2513: return(0);
2514: }
2516: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2517: {
2518: PetscInt i,nz = A->rmap->n*A->cmap->n;
2519: PetscScalar *aa;
2523: MatDenseGetArray(A,&aa);
2524: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2525: MatDenseRestoreArray(A,&aa);
2526: return(0);
2527: }
2529: /* ----------------------------------------------------------------*/
2530: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2531: {
2533: PetscInt m=A->rmap->n,n=B->cmap->n;
2534: PetscBool cisdense;
2537: MatSetSizes(C,m,n,m,n);
2538: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2539: if (!cisdense) {
2540: PetscBool flg;
2542: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2543: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2544: }
2545: MatSetUp(C);
2546: return(0);
2547: }
2549: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2550: {
2551: Mat_SeqDense *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data;
2552: PetscBLASInt m,n,k;
2553: const PetscScalar *av,*bv;
2554: PetscScalar *cv;
2555: PetscScalar _DOne=1.0,_DZero=0.0;
2556: PetscErrorCode ierr;
2559: PetscBLASIntCast(C->rmap->n,&m);
2560: PetscBLASIntCast(C->cmap->n,&n);
2561: PetscBLASIntCast(A->cmap->n,&k);
2562: if (!m || !n || !k) return(0);
2563: MatDenseGetArrayRead(A,&av);
2564: MatDenseGetArrayRead(B,&bv);
2565: MatDenseGetArrayWrite(C,&cv);
2566: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2567: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2568: MatDenseRestoreArrayRead(A,&av);
2569: MatDenseRestoreArrayRead(B,&bv);
2570: MatDenseRestoreArrayWrite(C,&cv);
2571: return(0);
2572: }
2574: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2575: {
2577: PetscInt m=A->rmap->n,n=B->rmap->n;
2578: PetscBool cisdense;
2581: MatSetSizes(C,m,n,m,n);
2582: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2583: if (!cisdense) {
2584: PetscBool flg;
2586: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2587: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2588: }
2589: MatSetUp(C);
2590: return(0);
2591: }
2593: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2594: {
2595: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2596: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2597: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2598: const PetscScalar *av,*bv;
2599: PetscScalar *cv;
2600: PetscBLASInt m,n,k;
2601: PetscScalar _DOne=1.0,_DZero=0.0;
2602: PetscErrorCode ierr;
2605: PetscBLASIntCast(C->rmap->n,&m);
2606: PetscBLASIntCast(C->cmap->n,&n);
2607: PetscBLASIntCast(A->cmap->n,&k);
2608: if (!m || !n || !k) return(0);
2609: MatDenseGetArrayRead(A,&av);
2610: MatDenseGetArrayRead(B,&bv);
2611: MatDenseGetArrayWrite(C,&cv);
2612: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2613: MatDenseRestoreArrayRead(A,&av);
2614: MatDenseRestoreArrayRead(B,&bv);
2615: MatDenseRestoreArrayWrite(C,&cv);
2616: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2617: return(0);
2618: }
2620: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2621: {
2623: PetscInt m=A->cmap->n,n=B->cmap->n;
2624: PetscBool cisdense;
2627: MatSetSizes(C,m,n,m,n);
2628: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2629: if (!cisdense) {
2630: PetscBool flg;
2632: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2633: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2634: }
2635: MatSetUp(C);
2636: return(0);
2637: }
2639: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2640: {
2641: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2642: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2643: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2644: const PetscScalar *av,*bv;
2645: PetscScalar *cv;
2646: PetscBLASInt m,n,k;
2647: PetscScalar _DOne=1.0,_DZero=0.0;
2648: PetscErrorCode ierr;
2651: PetscBLASIntCast(C->rmap->n,&m);
2652: PetscBLASIntCast(C->cmap->n,&n);
2653: PetscBLASIntCast(A->rmap->n,&k);
2654: if (!m || !n || !k) return(0);
2655: MatDenseGetArrayRead(A,&av);
2656: MatDenseGetArrayRead(B,&bv);
2657: MatDenseGetArrayWrite(C,&cv);
2658: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2659: MatDenseRestoreArrayRead(A,&av);
2660: MatDenseRestoreArrayRead(B,&bv);
2661: MatDenseRestoreArrayWrite(C,&cv);
2662: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2663: return(0);
2664: }
2666: /* ----------------------------------------------- */
2667: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2668: {
2670: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2671: C->ops->productsymbolic = MatProductSymbolic_AB;
2672: return(0);
2673: }
2675: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2676: {
2678: C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2679: C->ops->productsymbolic = MatProductSymbolic_AtB;
2680: return(0);
2681: }
2683: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2684: {
2686: C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2687: C->ops->productsymbolic = MatProductSymbolic_ABt;
2688: return(0);
2689: }
2691: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2692: {
2694: Mat_Product *product = C->product;
2697: switch (product->type) {
2698: case MATPRODUCT_AB:
2699: MatProductSetFromOptions_SeqDense_AB(C);
2700: break;
2701: case MATPRODUCT_AtB:
2702: MatProductSetFromOptions_SeqDense_AtB(C);
2703: break;
2704: case MATPRODUCT_ABt:
2705: MatProductSetFromOptions_SeqDense_ABt(C);
2706: break;
2707: default:
2708: break;
2709: }
2710: return(0);
2711: }
2712: /* ----------------------------------------------- */
2714: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2715: {
2716: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2717: PetscErrorCode ierr;
2718: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2719: PetscScalar *x;
2720: const PetscScalar *aa;
2723: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2724: VecGetArray(v,&x);
2725: VecGetLocalSize(v,&p);
2726: MatDenseGetArrayRead(A,&aa);
2727: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2728: for (i=0; i<m; i++) {
2729: x[i] = aa[i]; if (idx) idx[i] = 0;
2730: for (j=1; j<n; j++) {
2731: if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2732: }
2733: }
2734: MatDenseRestoreArrayRead(A,&aa);
2735: VecRestoreArray(v,&x);
2736: return(0);
2737: }
2739: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2740: {
2741: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2742: PetscErrorCode ierr;
2743: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2744: PetscScalar *x;
2745: PetscReal atmp;
2746: const PetscScalar *aa;
2749: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2750: VecGetArray(v,&x);
2751: VecGetLocalSize(v,&p);
2752: MatDenseGetArrayRead(A,&aa);
2753: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2754: for (i=0; i<m; i++) {
2755: x[i] = PetscAbsScalar(aa[i]);
2756: for (j=1; j<n; j++) {
2757: atmp = PetscAbsScalar(aa[i+a->lda*j]);
2758: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2759: }
2760: }
2761: MatDenseRestoreArrayRead(A,&aa);
2762: VecRestoreArray(v,&x);
2763: return(0);
2764: }
2766: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2767: {
2768: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2769: PetscErrorCode ierr;
2770: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2771: PetscScalar *x;
2772: const PetscScalar *aa;
2775: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2776: MatDenseGetArrayRead(A,&aa);
2777: VecGetArray(v,&x);
2778: VecGetLocalSize(v,&p);
2779: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2780: for (i=0; i<m; i++) {
2781: x[i] = aa[i]; if (idx) idx[i] = 0;
2782: for (j=1; j<n; j++) {
2783: if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2784: }
2785: }
2786: VecRestoreArray(v,&x);
2787: MatDenseRestoreArrayRead(A,&aa);
2788: return(0);
2789: }
2791: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2792: {
2793: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2794: PetscErrorCode ierr;
2795: PetscScalar *x;
2796: const PetscScalar *aa;
2799: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2800: MatDenseGetArrayRead(A,&aa);
2801: VecGetArray(v,&x);
2802: PetscArraycpy(x,aa+col*a->lda,A->rmap->n);
2803: VecRestoreArray(v,&x);
2804: MatDenseRestoreArrayRead(A,&aa);
2805: return(0);
2806: }
2808: PETSC_INTERN PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2809: {
2810: PetscErrorCode ierr;
2811: PetscInt i,j,m,n;
2812: const PetscScalar *a;
2815: MatGetSize(A,&m,&n);
2816: PetscArrayzero(norms,n);
2817: MatDenseGetArrayRead(A,&a);
2818: if (type == NORM_2) {
2819: for (i=0; i<n; i++) {
2820: for (j=0; j<m; j++) {
2821: norms[i] += PetscAbsScalar(a[j]*a[j]);
2822: }
2823: a += m;
2824: }
2825: } else if (type == NORM_1) {
2826: for (i=0; i<n; i++) {
2827: for (j=0; j<m; j++) {
2828: norms[i] += PetscAbsScalar(a[j]);
2829: }
2830: a += m;
2831: }
2832: } else if (type == NORM_INFINITY) {
2833: for (i=0; i<n; i++) {
2834: for (j=0; j<m; j++) {
2835: norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2836: }
2837: a += m;
2838: }
2839: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2840: MatDenseRestoreArrayRead(A,&a);
2841: if (type == NORM_2) {
2842: for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2843: }
2844: return(0);
2845: }
2847: static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2848: {
2850: PetscScalar *a;
2851: PetscInt lda,m,n,i,j;
2854: MatGetSize(x,&m,&n);
2855: MatDenseGetLDA(x,&lda);
2856: MatDenseGetArray(x,&a);
2857: for (j=0; j<n; j++) {
2858: for (i=0; i<m; i++) {
2859: PetscRandomGetValue(rctx,a+j*lda+i);
2860: }
2861: }
2862: MatDenseRestoreArray(x,&a);
2863: return(0);
2864: }
2866: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d)
2867: {
2869: *missing = PETSC_FALSE;
2870: return(0);
2871: }
2873: /* vals is not const */
2874: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2875: {
2877: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2878: PetscScalar *v;
2881: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2882: MatDenseGetArray(A,&v);
2883: *vals = v+col*a->lda;
2884: MatDenseRestoreArray(A,&v);
2885: return(0);
2886: }
2888: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2889: {
2891: *vals = NULL; /* user cannot accidently use the array later */
2892: return(0);
2893: }
2895: /* -------------------------------------------------------------------*/
2896: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2897: MatGetRow_SeqDense,
2898: MatRestoreRow_SeqDense,
2899: MatMult_SeqDense,
2900: /* 4*/ MatMultAdd_SeqDense,
2901: MatMultTranspose_SeqDense,
2902: MatMultTransposeAdd_SeqDense,
2903: NULL,
2904: NULL,
2905: NULL,
2906: /* 10*/ NULL,
2907: MatLUFactor_SeqDense,
2908: MatCholeskyFactor_SeqDense,
2909: MatSOR_SeqDense,
2910: MatTranspose_SeqDense,
2911: /* 15*/ MatGetInfo_SeqDense,
2912: MatEqual_SeqDense,
2913: MatGetDiagonal_SeqDense,
2914: MatDiagonalScale_SeqDense,
2915: MatNorm_SeqDense,
2916: /* 20*/ MatAssemblyBegin_SeqDense,
2917: MatAssemblyEnd_SeqDense,
2918: MatSetOption_SeqDense,
2919: MatZeroEntries_SeqDense,
2920: /* 24*/ MatZeroRows_SeqDense,
2921: NULL,
2922: NULL,
2923: NULL,
2924: NULL,
2925: /* 29*/ MatSetUp_SeqDense,
2926: NULL,
2927: NULL,
2928: NULL,
2929: NULL,
2930: /* 34*/ MatDuplicate_SeqDense,
2931: NULL,
2932: NULL,
2933: NULL,
2934: NULL,
2935: /* 39*/ MatAXPY_SeqDense,
2936: MatCreateSubMatrices_SeqDense,
2937: NULL,
2938: MatGetValues_SeqDense,
2939: MatCopy_SeqDense,
2940: /* 44*/ MatGetRowMax_SeqDense,
2941: MatScale_SeqDense,
2942: MatShift_Basic,
2943: NULL,
2944: MatZeroRowsColumns_SeqDense,
2945: /* 49*/ MatSetRandom_SeqDense,
2946: NULL,
2947: NULL,
2948: NULL,
2949: NULL,
2950: /* 54*/ NULL,
2951: NULL,
2952: NULL,
2953: NULL,
2954: NULL,
2955: /* 59*/ MatCreateSubMatrix_SeqDense,
2956: MatDestroy_SeqDense,
2957: MatView_SeqDense,
2958: NULL,
2959: NULL,
2960: /* 64*/ NULL,
2961: NULL,
2962: NULL,
2963: NULL,
2964: NULL,
2965: /* 69*/ MatGetRowMaxAbs_SeqDense,
2966: NULL,
2967: NULL,
2968: NULL,
2969: NULL,
2970: /* 74*/ NULL,
2971: NULL,
2972: NULL,
2973: NULL,
2974: NULL,
2975: /* 79*/ NULL,
2976: NULL,
2977: NULL,
2978: NULL,
2979: /* 83*/ MatLoad_SeqDense,
2980: MatIsSymmetric_SeqDense,
2981: MatIsHermitian_SeqDense,
2982: NULL,
2983: NULL,
2984: NULL,
2985: /* 89*/ NULL,
2986: NULL,
2987: MatMatMultNumeric_SeqDense_SeqDense,
2988: NULL,
2989: NULL,
2990: /* 94*/ NULL,
2991: NULL,
2992: NULL,
2993: MatMatTransposeMultNumeric_SeqDense_SeqDense,
2994: NULL,
2995: /* 99*/ MatProductSetFromOptions_SeqDense,
2996: NULL,
2997: NULL,
2998: MatConjugate_SeqDense,
2999: NULL,
3000: /*104*/ NULL,
3001: MatRealPart_SeqDense,
3002: MatImaginaryPart_SeqDense,
3003: NULL,
3004: NULL,
3005: /*109*/ NULL,
3006: NULL,
3007: MatGetRowMin_SeqDense,
3008: MatGetColumnVector_SeqDense,
3009: MatMissingDiagonal_SeqDense,
3010: /*114*/ NULL,
3011: NULL,
3012: NULL,
3013: NULL,
3014: NULL,
3015: /*119*/ NULL,
3016: NULL,
3017: NULL,
3018: NULL,
3019: NULL,
3020: /*124*/ NULL,
3021: MatGetColumnNorms_SeqDense,
3022: NULL,
3023: NULL,
3024: NULL,
3025: /*129*/ NULL,
3026: NULL,
3027: NULL,
3028: MatTransposeMatMultNumeric_SeqDense_SeqDense,
3029: NULL,
3030: /*134*/ NULL,
3031: NULL,
3032: NULL,
3033: NULL,
3034: NULL,
3035: /*139*/ NULL,
3036: NULL,
3037: NULL,
3038: NULL,
3039: NULL,
3040: MatCreateMPIMatConcatenateSeqMat_SeqDense,
3041: /*145*/ NULL,
3042: NULL,
3043: NULL
3044: };
3046: /*@C
3047: MatCreateSeqDense - Creates a sequential dense matrix that
3048: is stored in column major order (the usual Fortran 77 manner). Many
3049: of the matrix operations use the BLAS and LAPACK routines.
3051: Collective
3053: Input Parameters:
3054: + comm - MPI communicator, set to PETSC_COMM_SELF
3055: . m - number of rows
3056: . n - number of columns
3057: - data - optional location of matrix data in column major order. Set data=NULL for PETSc
3058: to control all matrix memory allocation.
3060: Output Parameter:
3061: . A - the matrix
3063: Notes:
3064: The data input variable is intended primarily for Fortran programmers
3065: who wish to allocate their own matrix memory space. Most users should
3066: set data=NULL.
3068: Level: intermediate
3070: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
3071: @*/
3072: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
3073: {
3077: MatCreate(comm,A);
3078: MatSetSizes(*A,m,n,m,n);
3079: MatSetType(*A,MATSEQDENSE);
3080: MatSeqDenseSetPreallocation(*A,data);
3081: return(0);
3082: }
3084: /*@C
3085: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
3087: Collective
3089: Input Parameters:
3090: + B - the matrix
3091: - data - the array (or NULL)
3093: Notes:
3094: The data input variable is intended primarily for Fortran programmers
3095: who wish to allocate their own matrix memory space. Most users should
3096: need not call this routine.
3098: Level: intermediate
3100: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatDenseSetLDA()
3102: @*/
3103: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
3104: {
3109: PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
3110: return(0);
3111: }
3113: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
3114: {
3115: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
3119: if (b->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3120: B->preallocated = PETSC_TRUE;
3122: PetscLayoutSetUp(B->rmap);
3123: PetscLayoutSetUp(B->cmap);
3125: if (b->lda <= 0) b->lda = B->rmap->n;
3127: PetscIntMultError(b->lda,B->cmap->n,NULL);
3128: if (!data) { /* petsc-allocated storage */
3129: if (!b->user_alloc) { PetscFree(b->v); }
3130: PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v);
3131: PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar));
3133: b->user_alloc = PETSC_FALSE;
3134: } else { /* user-allocated storage */
3135: if (!b->user_alloc) { PetscFree(b->v); }
3136: b->v = data;
3137: b->user_alloc = PETSC_TRUE;
3138: }
3139: B->assembled = PETSC_TRUE;
3140: return(0);
3141: }
3143: #if defined(PETSC_HAVE_ELEMENTAL)
3144: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
3145: {
3146: Mat mat_elemental;
3147: PetscErrorCode ierr;
3148: const PetscScalar *array;
3149: PetscScalar *v_colwise;
3150: PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;
3153: PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
3154: MatDenseGetArrayRead(A,&array);
3155: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
3156: k = 0;
3157: for (j=0; j<N; j++) {
3158: cols[j] = j;
3159: for (i=0; i<M; i++) {
3160: v_colwise[j*M+i] = array[k++];
3161: }
3162: }
3163: for (i=0; i<M; i++) {
3164: rows[i] = i;
3165: }
3166: MatDenseRestoreArrayRead(A,&array);
3168: MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
3169: MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
3170: MatSetType(mat_elemental,MATELEMENTAL);
3171: MatSetUp(mat_elemental);
3173: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
3174: MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
3175: MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
3176: MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
3177: PetscFree3(v_colwise,rows,cols);
3179: if (reuse == MAT_INPLACE_MATRIX) {
3180: MatHeaderReplace(A,&mat_elemental);
3181: } else {
3182: *newmat = mat_elemental;
3183: }
3184: return(0);
3185: }
3186: #endif
3188: PetscErrorCode MatDenseSetLDA_SeqDense(Mat B,PetscInt lda)
3189: {
3190: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
3191: PetscBool data;
3194: data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
3195: if (!b->user_alloc && data && b->lda!=lda) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"LDA cannot be changed after allocation of internal storage");
3196: 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);
3197: b->lda = lda;
3198: return(0);
3199: }
3201: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
3202: {
3204: PetscMPIInt size;
3207: MPI_Comm_size(comm,&size);
3208: if (size == 1) {
3209: if (scall == MAT_INITIAL_MATRIX) {
3210: MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
3211: } else {
3212: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
3213: }
3214: } else {
3215: MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
3216: }
3217: return(0);
3218: }
3220: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
3221: {
3222: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3226: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
3227: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3228: if (!a->cvec) {
3229: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3230: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3231: }
3232: a->vecinuse = col + 1;
3233: MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);
3234: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3235: *v = a->cvec;
3236: return(0);
3237: }
3239: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
3240: {
3241: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3245: if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
3246: if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
3247: a->vecinuse = 0;
3248: MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);
3249: VecResetArray(a->cvec);
3250: *v = NULL;
3251: return(0);
3252: }
3254: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
3255: {
3256: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3260: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
3261: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3262: if (!a->cvec) {
3263: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3264: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3265: }
3266: a->vecinuse = col + 1;
3267: MatDenseGetArrayRead(A,&a->ptrinuse);
3268: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3269: VecLockReadPush(a->cvec);
3270: *v = a->cvec;
3271: return(0);
3272: }
3274: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
3275: {
3276: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3280: if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
3281: if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
3282: a->vecinuse = 0;
3283: MatDenseRestoreArrayRead(A,&a->ptrinuse);
3284: VecLockReadPop(a->cvec);
3285: VecResetArray(a->cvec);
3286: *v = NULL;
3287: return(0);
3288: }
3290: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
3291: {
3292: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3296: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
3297: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3298: if (!a->cvec) {
3299: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3300: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3301: }
3302: a->vecinuse = col + 1;
3303: MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);
3304: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3305: *v = a->cvec;
3306: return(0);
3307: }
3309: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
3310: {
3311: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3315: if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
3316: if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
3317: a->vecinuse = 0;
3318: MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);
3319: VecResetArray(a->cvec);
3320: *v = NULL;
3321: return(0);
3322: }
3324: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3325: {
3326: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3330: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
3331: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3332: if (a->cmat && cend-cbegin != a->cmat->cmap->N) {
3333: MatDestroy(&a->cmat);
3334: }
3335: if (!a->cmat) {
3336: MatCreateDense(PetscObjectComm((PetscObject)A),A->rmap->n,PETSC_DECIDE,A->rmap->N,cend-cbegin,(PetscScalar*)a->v + (size_t)cbegin * (size_t)a->lda,&a->cmat);
3337: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat);
3338: } else {
3339: MatDensePlaceArray(a->cmat,a->v + (size_t)cbegin * (size_t)a->lda);
3340: }
3341: MatDenseSetLDA(a->cmat,a->lda);
3342: a->matinuse = cbegin + 1;
3343: *v = a->cmat;
3344: return(0);
3345: }
3347: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v)
3348: {
3349: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3353: if (!a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetSubMatrix() first");
3354: if (!a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column matrix");
3355: if (*v != a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not the matrix obtained from MatDenseGetSubMatrix()");
3356: a->matinuse = 0;
3357: MatDenseResetArray(a->cmat);
3358: *v = NULL;
3359: return(0);
3360: }
3362: /*MC
3363: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
3365: Options Database Keys:
3366: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
3368: Level: beginner
3370: .seealso: MatCreateSeqDense()
3372: M*/
3373: PetscErrorCode MatCreate_SeqDense(Mat B)
3374: {
3375: Mat_SeqDense *b;
3377: PetscMPIInt size;
3380: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
3381: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
3383: PetscNewLog(B,&b);
3384: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
3385: B->data = (void*)b;
3387: b->roworiented = PETSC_TRUE;
3389: PetscObjectComposeFunction((PetscObject)B,"MatQRFactor_C",MatQRFactor_SeqDense);
3390: PetscObjectComposeFunction((PetscObject)B,"MatQRFactorNumeric_C",MatQRFactorNumeric_SeqDense);
3391: PetscObjectComposeFunction((PetscObject)B,"MatQRFactorSymbolic_C",MatQRFactorSymbolic_SeqDense);
3392: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);
3393: PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense);
3394: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
3395: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
3396: PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
3397: PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
3398: PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense);
3399: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
3400: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);
3401: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense);
3402: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense);
3403: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
3404: #if defined(PETSC_HAVE_ELEMENTAL)
3405: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
3406: #endif
3407: #if defined(PETSC_HAVE_SCALAPACK)
3408: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_scalapack_C",MatConvert_Dense_ScaLAPACK);
3409: #endif
3410: #if defined(PETSC_HAVE_CUDA)
3411: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);
3412: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3413: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);
3414: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3415: #endif
3416: PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
3417: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);
3418: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);
3419: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
3420: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
3422: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
3423: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
3424: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);
3425: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);
3426: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);
3427: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);
3428: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);
3429: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);
3430: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense);
3431: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense);
3432: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
3433: return(0);
3434: }
3436: /*@C
3437: 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.
3439: Not Collective
3441: Input Parameters:
3442: + mat - a MATSEQDENSE or MATMPIDENSE matrix
3443: - col - column index
3445: Output Parameter:
3446: . vals - pointer to the data
3448: Level: intermediate
3450: .seealso: MatDenseRestoreColumn()
3451: @*/
3452: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
3453: {
3460: PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
3461: return(0);
3462: }
3464: /*@C
3465: MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn().
3467: Not Collective
3469: Input Parameter:
3470: . mat - a MATSEQDENSE or MATMPIDENSE matrix
3472: Output Parameter:
3473: . vals - pointer to the data
3475: Level: intermediate
3477: .seealso: MatDenseGetColumn()
3478: @*/
3479: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
3480: {
3486: PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
3487: return(0);
3488: }
3490: /*@C
3491: MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec.
3493: Collective
3495: Input Parameters:
3496: + mat - the Mat object
3497: - col - the column index
3499: Output Parameter:
3500: . v - the vector
3502: Notes:
3503: The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVec() when the vector is no longer needed.
3504: Use MatDenseGetColumnVecRead() to obtain read-only access or MatDenseGetColumnVecWrite() for write-only access.
3506: Level: intermediate
3508: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3509: @*/
3510: PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v)
3511: {
3519: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3520: if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3521: PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3522: return(0);
3523: }
3525: /*@C
3526: MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec().
3528: Collective
3530: Input Parameters:
3531: + mat - the Mat object
3532: . col - the column index
3533: - v - the Vec object
3535: Level: intermediate
3537: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3538: @*/
3539: PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v)
3540: {
3548: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3549: if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3550: PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3551: return(0);
3552: }
3554: /*@C
3555: MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec.
3557: Collective
3559: Input Parameters:
3560: + mat - the Mat object
3561: - col - the column index
3563: Output Parameter:
3564: . v - the vector
3566: Notes:
3567: The vector is owned by PETSc and users cannot modify it.
3568: Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed.
3569: Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access.
3571: Level: intermediate
3573: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3574: @*/
3575: PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v)
3576: {
3584: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3585: if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3586: PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3587: return(0);
3588: }
3590: /*@C
3591: MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead().
3593: Collective
3595: Input Parameters:
3596: + mat - the Mat object
3597: . col - the column index
3598: - v - the Vec object
3600: Level: intermediate
3602: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite()
3603: @*/
3604: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v)
3605: {
3613: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3614: if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3615: PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3616: return(0);
3617: }
3619: /*@C
3620: MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec.
3622: Collective
3624: Input Parameters:
3625: + mat - the Mat object
3626: - col - the column index
3628: Output Parameter:
3629: . v - the vector
3631: Notes:
3632: The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed.
3633: Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access.
3635: Level: intermediate
3637: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3638: @*/
3639: PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v)
3640: {
3648: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3649: if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3650: PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3651: return(0);
3652: }
3654: /*@C
3655: MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite().
3657: Collective
3659: Input Parameters:
3660: + mat - the Mat object
3661: . col - the column index
3662: - v - the Vec object
3664: Level: intermediate
3666: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead()
3667: @*/
3668: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v)
3669: {
3677: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3678: if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N);
3679: PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3680: return(0);
3681: }
3683: /*@C
3684: MatDenseGetSubMatrix - Gives access to a block of columns of a dense matrix, represented as a Mat.
3686: Collective
3688: Input Parameters:
3689: + mat - the Mat object
3690: . cbegin - the first index in the block
3691: - cend - the last index in the block
3693: Output Parameter:
3694: . v - the matrix
3696: Notes:
3697: The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed.
3699: Level: intermediate
3701: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseRestoreSubMatrix()
3702: @*/
3703: PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3704: {
3713: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3714: if (cbegin < 0 || cbegin > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cbegin %D, should be in [0,%D)",cbegin,A->cmap->N);
3715: if (cend < cbegin || cend > A->cmap->N) SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cend %D, should be in [%D,%D)",cend,cbegin,A->cmap->N);
3716: PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));
3717: return(0);
3718: }
3720: /*@C
3721: MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix().
3723: Collective
3725: Input Parameters:
3726: + mat - the Mat object
3727: - v - the Mat object
3729: Level: intermediate
3731: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseGetSubMatrix()
3732: @*/
3733: PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v)
3734: {
3741: PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v));
3742: return(0);
3743: }