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 ldx, 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 ldx, 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 ldx, 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,&ldx,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,&ldx,&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*ldx + 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 ldx, 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,&ldx,&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,&ldx,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 *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
670: {
671: PetscErrorCode ierr;
672: const PetscScalar *b;
673: PetscScalar *y;
674: PetscInt n, _ldb, _ldx;
675: PetscBLASInt nrhs=0,m=0,k=0,ldb=0,ldx=0,ldy=0;
678: *_ldy=0; *_m=0; *_nrhs=0; *_k=0; *_y = NULL;
679: PetscBLASIntCast(A->rmap->n,&m);
680: PetscBLASIntCast(A->cmap->n,&k);
681: MatGetSize(B,NULL,&n);
682: PetscBLASIntCast(n,&nrhs);
683: MatDenseGetLDA(B,&_ldb);
684: PetscBLASIntCast(_ldb, &ldb);
685: MatDenseGetLDA(X,&_ldx);
686: PetscBLASIntCast(_ldx, &ldx);
687: if (ldx < m) {
688: MatDenseGetArrayRead(B,&b);
689: PetscMalloc1(nrhs * m, &y);
690: if (ldb == m) {
691: PetscArraycpy(y,b,ldb*nrhs);
692: } else {
693: for (PetscInt j = 0; j < nrhs; j++) {
694: PetscArraycpy(&y[j*m],&b[j*ldb],m);
695: }
696: }
697: ldy = m;
698: MatDenseRestoreArrayRead(B,&b);
699: } else {
700: if (ldb == ldx) {
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*ldx],&b[j*ldb],m);
708: }
709: MatDenseRestoreArrayRead(B,&b);
710: }
711: ldy = ldx;
712: }
713: *_y = y;
714: *_ldy = ldy;
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 *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
722: {
723: PetscScalar *y;
724: PetscInt _ldx;
725: PetscBLASInt k,ldy,nrhs,ldx=0;
726: PetscErrorCode ierr;
729: y = *_y;
730: *_y = NULL;
731: k = *_k;
732: ldy = *_ldy;
733: nrhs = *_nrhs;
734: MatDenseGetLDA(X,&_ldx);
735: PetscBLASIntCast(_ldx, &ldx);
736: if (ldx != ldy) {
737: PetscScalar *xv;
738: MatDenseGetArray(X,&xv);
739: for (PetscInt j = 0; j < nrhs; j++) {
740: PetscArraycpy(&xv[j*ldx],&y[j*ldy],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, ldy, nrhs;
757: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
758: MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_FALSE);
759: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &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, ldy, nrhs;
770: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
771: MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_TRUE);
772: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &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, ldy, nrhs;
783: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
784: MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_FALSE);
785: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &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, ldy, nrhs;
796: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
797: MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_TRUE);
798: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &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, ldy, nrhs;
809: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
810: MatSolve_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k);
811: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &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, ldy, nrhs;
822: MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k);
823: MatSolveTranspose_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k);
824: MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &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: return(0);
885: }
887: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
888: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
889: {
890: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
892: PetscBLASInt info,n;
895: PetscBLASIntCast(A->cmap->n,&n);
896: if (!A->rmap->n || !A->cmap->n) return(0);
897: if (A->spd) {
898: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
899: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
900: PetscFPTrapPop();
901: #if defined(PETSC_USE_COMPLEX)
902: } else if (A->hermitian) {
903: if (!mat->pivots) {
904: PetscMalloc1(A->rmap->n,&mat->pivots);
905: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
906: }
907: if (!mat->fwork) {
908: PetscScalar dummy;
910: mat->lfwork = -1;
911: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
912: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
913: PetscFPTrapPop();
914: mat->lfwork = (PetscInt)PetscRealPart(dummy);
915: PetscMalloc1(mat->lfwork,&mat->fwork);
916: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
917: }
918: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
919: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
920: PetscFPTrapPop();
921: #endif
922: } else { /* symmetric case */
923: if (!mat->pivots) {
924: PetscMalloc1(A->rmap->n,&mat->pivots);
925: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
926: }
927: if (!mat->fwork) {
928: PetscScalar dummy;
930: mat->lfwork = -1;
931: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
932: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
933: PetscFPTrapPop();
934: mat->lfwork = (PetscInt)PetscRealPart(dummy);
935: PetscMalloc1(mat->lfwork,&mat->fwork);
936: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
937: }
938: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
939: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
940: PetscFPTrapPop();
941: }
942: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
944: A->ops->solve = MatSolve_SeqDense_Cholesky;
945: A->ops->matsolve = MatMatSolve_SeqDense_Cholesky;
946: A->ops->solvetranspose = MatSolveTranspose_SeqDense_Cholesky;
947: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky;
948: A->factortype = MAT_FACTOR_CHOLESKY;
950: PetscFree(A->solvertype);
951: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
953: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
954: return(0);
955: }
957: static PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
958: {
960: MatFactorInfo info;
963: info.fill = 1.0;
965: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
966: (*fact->ops->choleskyfactor)(fact,NULL,&info);
967: return(0);
968: }
970: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
971: {
973: fact->assembled = PETSC_TRUE;
974: fact->preallocated = PETSC_TRUE;
975: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
976: return(0);
977: }
979: PetscErrorCode MatQRFactor_SeqDense(Mat A,IS col,const MatFactorInfo *minfo)
980: {
981: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
983: PetscBLASInt n,m,info, min, max;
986: PetscBLASIntCast(A->cmap->n,&n);
987: PetscBLASIntCast(A->rmap->n,&m);
988: max = PetscMax(m, n);
989: min = PetscMin(m, n);
990: if (!mat->tau) {
991: PetscMalloc1(min,&mat->tau);
992: PetscLogObjectMemory((PetscObject)A,min*sizeof(PetscScalar));
993: }
994: if (!mat->pivots) {
995: PetscMalloc1(n,&mat->pivots);
996: PetscLogObjectMemory((PetscObject)A,n*sizeof(PetscScalar));
997: }
998: if (!mat->qrrhs) {
999: MatCreateVecs(A, NULL, &(mat->qrrhs));
1000: }
1001: if (!A->rmap->n || !A->cmap->n) return(0);
1002: if (!mat->fwork) {
1003: PetscScalar dummy;
1005: mat->lfwork = -1;
1006: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1007: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,&dummy,&mat->lfwork,&info));
1008: PetscFPTrapPop();
1009: mat->lfwork = (PetscInt)PetscRealPart(dummy);
1010: PetscMalloc1(mat->lfwork,&mat->fwork);
1011: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
1012: }
1013: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1014: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,mat->v,&mat->lda,mat->tau,mat->fwork,&mat->lfwork,&info));
1015: PetscFPTrapPop();
1016: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to QR factorization");
1017: // 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
1018: mat->rank = min;
1020: A->ops->solve = MatSolve_SeqDense_QR;
1021: A->ops->matsolve = MatMatSolve_SeqDense_QR;
1022: A->factortype = MAT_FACTOR_QR;
1023: if (m == n) {
1024: A->ops->solvetranspose = MatSolveTranspose_SeqDense_QR;
1025: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
1026: }
1028: PetscFree(A->solvertype);
1029: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
1031: PetscLogFlops(2.0*min*min*(max-min/3.0));
1032: return(0);
1033: }
1035: static PetscErrorCode MatQRFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
1036: {
1038: MatFactorInfo info;
1041: info.fill = 1.0;
1043: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
1044: PetscUseMethod(fact,"MatQRFactor_C",(Mat,IS,const MatFactorInfo *),(fact,NULL,&info));
1045: return(0);
1046: }
1048: PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
1049: {
1053: fact->assembled = PETSC_TRUE;
1054: fact->preallocated = PETSC_TRUE;
1055: PetscObjectComposeFunction((PetscObject)fact,"MatQRFactorNumeric_C",MatQRFactorNumeric_SeqDense);
1056: return(0);
1057: }
1059: /* uses LAPACK */
1060: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
1061: {
1065: MatCreate(PetscObjectComm((PetscObject)A),fact);
1066: MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
1067: MatSetType(*fact,MATDENSE);
1068: (*fact)->trivialsymbolic = PETSC_TRUE;
1069: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
1070: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
1071: (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
1072: } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
1073: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
1074: } else if (ftype == MAT_FACTOR_QR) {
1075: PetscObjectComposeFunction((PetscObject)(*fact),"MatQRFactorSymbolic_C",MatQRFactorSymbolic_SeqDense);
1076: }
1077: (*fact)->factortype = ftype;
1079: PetscFree((*fact)->solvertype);
1080: PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
1081: PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_LU]);
1082: PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_ILU]);
1083: PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_CHOLESKY]);
1084: PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&(*fact)->preferredordering[MAT_FACTOR_ICC]);
1085: return(0);
1086: }
1088: /* ------------------------------------------------------------------*/
1089: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
1090: {
1091: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1092: PetscScalar *x,*v = mat->v,zero = 0.0,xt;
1093: const PetscScalar *b;
1094: PetscErrorCode ierr;
1095: PetscInt m = A->rmap->n,i;
1096: PetscBLASInt o = 1,bm = 0;
1099: #if defined(PETSC_HAVE_CUDA)
1100: if (A->offloadmask == PETSC_OFFLOAD_GPU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented");
1101: #endif
1102: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
1103: PetscBLASIntCast(m,&bm);
1104: if (flag & SOR_ZERO_INITIAL_GUESS) {
1105: /* this is a hack fix, should have another version without the second BLASdotu */
1106: VecSet(xx,zero);
1107: }
1108: VecGetArray(xx,&x);
1109: VecGetArrayRead(bb,&b);
1110: its = its*lits;
1111: if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
1112: while (its--) {
1113: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1114: for (i=0; i<m; i++) {
1115: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
1116: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
1117: }
1118: }
1119: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1120: for (i=m-1; i>=0; 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: }
1126: VecRestoreArrayRead(bb,&b);
1127: VecRestoreArray(xx,&x);
1128: return(0);
1129: }
1131: /* -----------------------------------------------------------------*/
1132: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
1133: {
1134: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1135: const PetscScalar *v = mat->v,*x;
1136: PetscScalar *y;
1137: PetscErrorCode ierr;
1138: PetscBLASInt m, n,_One=1;
1139: PetscScalar _DOne=1.0,_DZero=0.0;
1142: PetscBLASIntCast(A->rmap->n,&m);
1143: PetscBLASIntCast(A->cmap->n,&n);
1144: VecGetArrayRead(xx,&x);
1145: VecGetArrayWrite(yy,&y);
1146: if (!A->rmap->n || !A->cmap->n) {
1147: PetscBLASInt i;
1148: for (i=0; i<n; i++) y[i] = 0.0;
1149: } else {
1150: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
1151: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
1152: }
1153: VecRestoreArrayRead(xx,&x);
1154: VecRestoreArrayWrite(yy,&y);
1155: return(0);
1156: }
1158: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
1159: {
1160: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1161: PetscScalar *y,_DOne=1.0,_DZero=0.0;
1162: PetscErrorCode ierr;
1163: PetscBLASInt m, n, _One=1;
1164: const PetscScalar *v = mat->v,*x;
1167: PetscBLASIntCast(A->rmap->n,&m);
1168: PetscBLASIntCast(A->cmap->n,&n);
1169: VecGetArrayRead(xx,&x);
1170: VecGetArrayWrite(yy,&y);
1171: if (!A->rmap->n || !A->cmap->n) {
1172: PetscBLASInt i;
1173: for (i=0; i<m; i++) y[i] = 0.0;
1174: } else {
1175: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
1176: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
1177: }
1178: VecRestoreArrayRead(xx,&x);
1179: VecRestoreArrayWrite(yy,&y);
1180: return(0);
1181: }
1183: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
1184: {
1185: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1186: const PetscScalar *v = mat->v,*x;
1187: PetscScalar *y,_DOne=1.0;
1188: PetscErrorCode ierr;
1189: PetscBLASInt m, n, _One=1;
1192: PetscBLASIntCast(A->rmap->n,&m);
1193: PetscBLASIntCast(A->cmap->n,&n);
1194: VecCopy(zz,yy);
1195: if (!A->rmap->n || !A->cmap->n) return(0);
1196: VecGetArrayRead(xx,&x);
1197: VecGetArray(yy,&y);
1198: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
1199: VecRestoreArrayRead(xx,&x);
1200: VecRestoreArray(yy,&y);
1201: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
1202: return(0);
1203: }
1205: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
1206: {
1207: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1208: const PetscScalar *v = mat->v,*x;
1209: PetscScalar *y;
1210: PetscErrorCode ierr;
1211: PetscBLASInt m, n, _One=1;
1212: PetscScalar _DOne=1.0;
1215: PetscBLASIntCast(A->rmap->n,&m);
1216: PetscBLASIntCast(A->cmap->n,&n);
1217: VecCopy(zz,yy);
1218: if (!A->rmap->n || !A->cmap->n) return(0);
1219: VecGetArrayRead(xx,&x);
1220: VecGetArray(yy,&y);
1221: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
1222: VecRestoreArrayRead(xx,&x);
1223: VecRestoreArray(yy,&y);
1224: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
1225: return(0);
1226: }
1228: /* -----------------------------------------------------------------*/
1229: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
1230: {
1231: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1233: PetscInt i;
1236: *ncols = A->cmap->n;
1237: if (cols) {
1238: PetscMalloc1(A->cmap->n,cols);
1239: for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
1240: }
1241: if (vals) {
1242: const PetscScalar *v;
1244: MatDenseGetArrayRead(A,&v);
1245: PetscMalloc1(A->cmap->n,vals);
1246: v += row;
1247: for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
1248: MatDenseRestoreArrayRead(A,&v);
1249: }
1250: return(0);
1251: }
1253: static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
1254: {
1258: if (ncols) *ncols = 0;
1259: if (cols) {PetscFree(*cols);}
1260: if (vals) {PetscFree(*vals); }
1261: return(0);
1262: }
1263: /* ----------------------------------------------------------------*/
1264: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
1265: {
1266: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1267: PetscScalar *av;
1268: PetscInt i,j,idx=0;
1269: #if defined(PETSC_HAVE_CUDA)
1270: PetscOffloadMask oldf;
1271: #endif
1272: PetscErrorCode ierr;
1275: MatDenseGetArray(A,&av);
1276: if (!mat->roworiented) {
1277: if (addv == INSERT_VALUES) {
1278: for (j=0; j<n; j++) {
1279: if (indexn[j] < 0) {idx += m; continue;}
1280: 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);
1281: for (i=0; i<m; i++) {
1282: if (indexm[i] < 0) {idx++; continue;}
1283: 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);
1284: av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1285: }
1286: }
1287: } else {
1288: for (j=0; j<n; j++) {
1289: if (indexn[j] < 0) {idx += m; continue;}
1290: 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);
1291: for (i=0; i<m; i++) {
1292: if (indexm[i] < 0) {idx++; continue;}
1293: 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);
1294: av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1295: }
1296: }
1297: }
1298: } else {
1299: if (addv == INSERT_VALUES) {
1300: for (i=0; i<m; i++) {
1301: if (indexm[i] < 0) { idx += n; continue;}
1302: 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);
1303: for (j=0; j<n; j++) {
1304: if (indexn[j] < 0) { idx++; continue;}
1305: 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);
1306: av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
1307: }
1308: }
1309: } else {
1310: for (i=0; i<m; i++) {
1311: if (indexm[i] < 0) { idx += n; continue;}
1312: 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);
1313: for (j=0; j<n; j++) {
1314: if (indexn[j] < 0) { idx++; continue;}
1315: 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);
1316: av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
1317: }
1318: }
1319: }
1320: }
1321: /* hack to prevent unneeded copy to the GPU while returning the array */
1322: #if defined(PETSC_HAVE_CUDA)
1323: oldf = A->offloadmask;
1324: A->offloadmask = PETSC_OFFLOAD_GPU;
1325: #endif
1326: MatDenseRestoreArray(A,&av);
1327: #if defined(PETSC_HAVE_CUDA)
1328: A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1329: #endif
1330: return(0);
1331: }
1333: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1334: {
1335: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1336: const PetscScalar *vv;
1337: PetscInt i,j;
1338: PetscErrorCode ierr;
1341: MatDenseGetArrayRead(A,&vv);
1342: /* row-oriented output */
1343: for (i=0; i<m; i++) {
1344: if (indexm[i] < 0) {v += n;continue;}
1345: 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);
1346: for (j=0; j<n; j++) {
1347: if (indexn[j] < 0) {v++; continue;}
1348: 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);
1349: *v++ = vv[indexn[j]*mat->lda + indexm[i]];
1350: }
1351: }
1352: MatDenseRestoreArrayRead(A,&vv);
1353: return(0);
1354: }
1356: /* -----------------------------------------------------------------*/
1358: PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer)
1359: {
1360: PetscErrorCode ierr;
1361: PetscBool skipHeader;
1362: PetscViewerFormat format;
1363: PetscInt header[4],M,N,m,lda,i,j,k;
1364: const PetscScalar *v;
1365: PetscScalar *vwork;
1368: PetscViewerSetUp(viewer);
1369: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
1370: PetscViewerGetFormat(viewer,&format);
1371: if (skipHeader) format = PETSC_VIEWER_NATIVE;
1373: MatGetSize(mat,&M,&N);
1375: /* write matrix header */
1376: header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N;
1377: header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N;
1378: if (!skipHeader) {PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);}
1380: MatGetLocalSize(mat,&m,NULL);
1381: if (format != PETSC_VIEWER_NATIVE) {
1382: PetscInt nnz = m*N, *iwork;
1383: /* store row lengths for each row */
1384: PetscMalloc1(nnz,&iwork);
1385: for (i=0; i<m; i++) iwork[i] = N;
1386: PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1387: /* store column indices (zero start index) */
1388: for (k=0, i=0; i<m; i++)
1389: for (j=0; j<N; j++, k++)
1390: iwork[k] = j;
1391: PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1392: PetscFree(iwork);
1393: }
1394: /* store matrix values as a dense matrix in row major order */
1395: PetscMalloc1(m*N,&vwork);
1396: MatDenseGetArrayRead(mat,&v);
1397: MatDenseGetLDA(mat,&lda);
1398: for (k=0, i=0; i<m; i++)
1399: for (j=0; j<N; j++, k++)
1400: vwork[k] = v[i+lda*j];
1401: MatDenseRestoreArrayRead(mat,&v);
1402: PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1403: PetscFree(vwork);
1404: return(0);
1405: }
1407: PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer)
1408: {
1410: PetscBool skipHeader;
1411: PetscInt header[4],M,N,m,nz,lda,i,j,k;
1412: PetscInt rows,cols;
1413: PetscScalar *v,*vwork;
1416: PetscViewerSetUp(viewer);
1417: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
1419: if (!skipHeader) {
1420: PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
1421: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file");
1422: M = header[1]; N = header[2];
1423: if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M);
1424: if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N);
1425: nz = header[3];
1426: if (nz != MATRIX_BINARY_FORMAT_DENSE && nz < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Unknown matrix format %D in file",nz);
1427: } else {
1428: MatGetSize(mat,&M,&N);
1429: 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");
1430: nz = MATRIX_BINARY_FORMAT_DENSE;
1431: }
1433: /* setup global sizes if not set */
1434: if (mat->rmap->N < 0) mat->rmap->N = M;
1435: if (mat->cmap->N < 0) mat->cmap->N = N;
1436: MatSetUp(mat);
1437: /* check if global sizes are correct */
1438: MatGetSize(mat,&rows,&cols);
1439: 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);
1441: MatGetSize(mat,NULL,&N);
1442: MatGetLocalSize(mat,&m,NULL);
1443: MatDenseGetArray(mat,&v);
1444: MatDenseGetLDA(mat,&lda);
1445: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */
1446: PetscInt nnz = m*N;
1447: /* read in matrix values */
1448: PetscMalloc1(nnz,&vwork);
1449: PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1450: /* store values in column major order */
1451: for (j=0; j<N; j++)
1452: for (i=0; i<m; i++)
1453: v[i+lda*j] = vwork[i*N+j];
1454: PetscFree(vwork);
1455: } else { /* matrix in file is sparse format */
1456: PetscInt nnz = 0, *rlens, *icols;
1457: /* read in row lengths */
1458: PetscMalloc1(m,&rlens);
1459: PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1460: for (i=0; i<m; i++) nnz += rlens[i];
1461: /* read in column indices and values */
1462: PetscMalloc2(nnz,&icols,nnz,&vwork);
1463: PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1464: PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1465: /* store values in column major order */
1466: for (k=0, i=0; i<m; i++)
1467: for (j=0; j<rlens[i]; j++, k++)
1468: v[i+lda*icols[k]] = vwork[k];
1469: PetscFree(rlens);
1470: PetscFree2(icols,vwork);
1471: }
1472: MatDenseRestoreArray(mat,&v);
1473: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
1474: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
1475: return(0);
1476: }
1478: PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1479: {
1480: PetscBool isbinary, ishdf5;
1486: /* force binary viewer to load .info file if it has not yet done so */
1487: PetscViewerSetUp(viewer);
1488: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1489: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);
1490: if (isbinary) {
1491: MatLoad_Dense_Binary(newMat,viewer);
1492: } else if (ishdf5) {
1493: #if defined(PETSC_HAVE_HDF5)
1494: MatLoad_Dense_HDF5(newMat,viewer);
1495: #else
1496: SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1497: #endif
1498: } else {
1499: 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);
1500: }
1501: return(0);
1502: }
1504: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1505: {
1506: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1507: PetscErrorCode ierr;
1508: PetscInt i,j;
1509: const char *name;
1510: PetscScalar *v,*av;
1511: PetscViewerFormat format;
1512: #if defined(PETSC_USE_COMPLEX)
1513: PetscBool allreal = PETSC_TRUE;
1514: #endif
1517: MatDenseGetArrayRead(A,(const PetscScalar**)&av);
1518: PetscViewerGetFormat(viewer,&format);
1519: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1520: return(0); /* do nothing for now */
1521: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1522: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1523: for (i=0; i<A->rmap->n; i++) {
1524: v = av + i;
1525: PetscViewerASCIIPrintf(viewer,"row %D:",i);
1526: for (j=0; j<A->cmap->n; j++) {
1527: #if defined(PETSC_USE_COMPLEX)
1528: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1529: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1530: } else if (PetscRealPart(*v)) {
1531: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1532: }
1533: #else
1534: if (*v) {
1535: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1536: }
1537: #endif
1538: v += a->lda;
1539: }
1540: PetscViewerASCIIPrintf(viewer,"\n");
1541: }
1542: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1543: } else {
1544: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1545: #if defined(PETSC_USE_COMPLEX)
1546: /* determine if matrix has all real values */
1547: v = av;
1548: for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1549: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1550: }
1551: #endif
1552: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1553: PetscObjectGetName((PetscObject)A,&name);
1554: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1555: PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1556: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1557: }
1559: for (i=0; i<A->rmap->n; i++) {
1560: v = av + i;
1561: for (j=0; j<A->cmap->n; j++) {
1562: #if defined(PETSC_USE_COMPLEX)
1563: if (allreal) {
1564: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1565: } else {
1566: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1567: }
1568: #else
1569: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1570: #endif
1571: v += a->lda;
1572: }
1573: PetscViewerASCIIPrintf(viewer,"\n");
1574: }
1575: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1576: PetscViewerASCIIPrintf(viewer,"];\n");
1577: }
1578: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1579: }
1580: MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);
1581: PetscViewerFlush(viewer);
1582: return(0);
1583: }
1585: #include <petscdraw.h>
1586: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1587: {
1588: Mat A = (Mat) Aa;
1589: PetscErrorCode ierr;
1590: PetscInt m = A->rmap->n,n = A->cmap->n,i,j;
1591: int color = PETSC_DRAW_WHITE;
1592: const PetscScalar *v;
1593: PetscViewer viewer;
1594: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1595: PetscViewerFormat format;
1598: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1599: PetscViewerGetFormat(viewer,&format);
1600: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
1602: /* Loop over matrix elements drawing boxes */
1603: MatDenseGetArrayRead(A,&v);
1604: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1605: PetscDrawCollectiveBegin(draw);
1606: /* Blue for negative and Red for positive */
1607: for (j = 0; j < n; j++) {
1608: x_l = j; x_r = x_l + 1.0;
1609: for (i = 0; i < m; i++) {
1610: y_l = m - i - 1.0;
1611: y_r = y_l + 1.0;
1612: if (PetscRealPart(v[j*m+i]) > 0.) color = PETSC_DRAW_RED;
1613: else if (PetscRealPart(v[j*m+i]) < 0.) color = PETSC_DRAW_BLUE;
1614: else continue;
1615: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1616: }
1617: }
1618: PetscDrawCollectiveEnd(draw);
1619: } else {
1620: /* use contour shading to indicate magnitude of values */
1621: /* first determine max of all nonzero values */
1622: PetscReal minv = 0.0, maxv = 0.0;
1623: PetscDraw popup;
1625: for (i=0; i < m*n; i++) {
1626: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1627: }
1628: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1629: PetscDrawGetPopup(draw,&popup);
1630: PetscDrawScalePopup(popup,minv,maxv);
1632: PetscDrawCollectiveBegin(draw);
1633: for (j=0; j<n; j++) {
1634: x_l = j;
1635: x_r = x_l + 1.0;
1636: for (i=0; i<m; i++) {
1637: y_l = m - i - 1.0;
1638: y_r = y_l + 1.0;
1639: color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1640: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1641: }
1642: }
1643: PetscDrawCollectiveEnd(draw);
1644: }
1645: MatDenseRestoreArrayRead(A,&v);
1646: return(0);
1647: }
1649: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1650: {
1651: PetscDraw draw;
1652: PetscBool isnull;
1653: PetscReal xr,yr,xl,yl,h,w;
1657: PetscViewerDrawGetDraw(viewer,0,&draw);
1658: PetscDrawIsNull(draw,&isnull);
1659: if (isnull) return(0);
1661: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1662: xr += w; yr += h; xl = -w; yl = -h;
1663: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1664: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1665: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1666: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1667: PetscDrawSave(draw);
1668: return(0);
1669: }
1671: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1672: {
1674: PetscBool iascii,isbinary,isdraw;
1677: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1678: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1679: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1681: if (iascii) {
1682: MatView_SeqDense_ASCII(A,viewer);
1683: } else if (isbinary) {
1684: MatView_Dense_Binary(A,viewer);
1685: } else if (isdraw) {
1686: MatView_SeqDense_Draw(A,viewer);
1687: }
1688: return(0);
1689: }
1691: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar *array)
1692: {
1693: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1696: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1697: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1698: if (a->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreArray() first");
1699: a->unplacedarray = a->v;
1700: a->unplaced_user_alloc = a->user_alloc;
1701: a->v = (PetscScalar*) array;
1702: a->user_alloc = PETSC_TRUE;
1703: #if defined(PETSC_HAVE_CUDA)
1704: A->offloadmask = PETSC_OFFLOAD_CPU;
1705: #endif
1706: return(0);
1707: }
1709: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1710: {
1711: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1714: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1715: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1716: a->v = a->unplacedarray;
1717: a->user_alloc = a->unplaced_user_alloc;
1718: a->unplacedarray = NULL;
1719: #if defined(PETSC_HAVE_CUDA)
1720: A->offloadmask = PETSC_OFFLOAD_CPU;
1721: #endif
1722: return(0);
1723: }
1725: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A,const PetscScalar *array)
1726: {
1727: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1731: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1732: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1733: if (!a->user_alloc) { PetscFree(a->v); }
1734: a->v = (PetscScalar*) array;
1735: a->user_alloc = PETSC_FALSE;
1736: #if defined(PETSC_HAVE_CUDA)
1737: A->offloadmask = PETSC_OFFLOAD_CPU;
1738: #endif
1739: return(0);
1740: }
1742: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1743: {
1744: Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
1748: #if defined(PETSC_USE_LOG)
1749: PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1750: #endif
1751: VecDestroy(&(l->qrrhs));
1752: PetscFree(l->tau);
1753: PetscFree(l->pivots);
1754: PetscFree(l->fwork);
1755: MatDestroy(&l->ptapwork);
1756: if (!l->user_alloc) {PetscFree(l->v);}
1757: if (!l->unplaced_user_alloc) {PetscFree(l->unplacedarray);}
1758: if (l->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1759: if (l->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1760: VecDestroy(&l->cvec);
1761: MatDestroy(&l->cmat);
1762: PetscFree(mat->data);
1764: PetscObjectChangeTypeName((PetscObject)mat,NULL);
1765: PetscObjectComposeFunction((PetscObject)mat,"MatQRFactor_C",NULL);
1766: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);
1767: PetscObjectComposeFunction((PetscObject)mat,"MatDenseSetLDA_C",NULL);
1768: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1769: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1770: PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1771: PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1772: PetscObjectComposeFunction((PetscObject)mat,"MatDenseReplaceArray_C",NULL);
1773: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);
1774: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);
1775: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayWrite_C",NULL);
1776: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayWrite_C",NULL);
1777: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1778: #if defined(PETSC_HAVE_ELEMENTAL)
1779: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1780: #endif
1781: #if defined(PETSC_HAVE_SCALAPACK)
1782: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_scalapack_C",NULL);
1783: #endif
1784: #if defined(PETSC_HAVE_CUDA)
1785: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);
1786: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL);
1787: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL);
1788: #endif
1789: PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1790: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL);
1791: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL);
1792: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL);
1793: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL);
1795: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
1796: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
1797: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVec_C",NULL);
1798: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVec_C",NULL);
1799: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecRead_C",NULL);
1800: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecRead_C",NULL);
1801: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecWrite_C",NULL);
1802: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecWrite_C",NULL);
1803: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetSubMatrix_C",NULL);
1804: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreSubMatrix_C",NULL);
1805: return(0);
1806: }
1808: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1809: {
1810: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1812: PetscInt k,j,m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1813: PetscScalar *v,tmp;
1816: if (reuse == MAT_INPLACE_MATRIX) {
1817: if (m == n) { /* in place transpose */
1818: MatDenseGetArray(A,&v);
1819: for (j=0; j<m; j++) {
1820: for (k=0; k<j; k++) {
1821: tmp = v[j + k*M];
1822: v[j + k*M] = v[k + j*M];
1823: v[k + j*M] = tmp;
1824: }
1825: }
1826: MatDenseRestoreArray(A,&v);
1827: } else { /* reuse memory, temporary allocates new memory */
1828: PetscScalar *v2;
1829: PetscLayout tmplayout;
1831: PetscMalloc1((size_t)m*n,&v2);
1832: MatDenseGetArray(A,&v);
1833: for (j=0; j<n; j++) {
1834: for (k=0; k<m; k++) v2[j + (size_t)k*n] = v[k + (size_t)j*M];
1835: }
1836: PetscArraycpy(v,v2,(size_t)m*n);
1837: PetscFree(v2);
1838: MatDenseRestoreArray(A,&v);
1839: /* cleanup size dependent quantities */
1840: VecDestroy(&mat->cvec);
1841: MatDestroy(&mat->cmat);
1842: PetscFree(mat->pivots);
1843: PetscFree(mat->fwork);
1844: MatDestroy(&mat->ptapwork);
1845: /* swap row/col layouts */
1846: mat->lda = n;
1847: tmplayout = A->rmap;
1848: A->rmap = A->cmap;
1849: A->cmap = tmplayout;
1850: }
1851: } else { /* out-of-place transpose */
1852: Mat tmat;
1853: Mat_SeqDense *tmatd;
1854: PetscScalar *v2;
1855: PetscInt M2;
1857: if (reuse == MAT_INITIAL_MATRIX) {
1858: MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1859: MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1860: MatSetType(tmat,((PetscObject)A)->type_name);
1861: MatSeqDenseSetPreallocation(tmat,NULL);
1862: } else tmat = *matout;
1864: MatDenseGetArrayRead(A,(const PetscScalar**)&v);
1865: MatDenseGetArray(tmat,&v2);
1866: tmatd = (Mat_SeqDense*)tmat->data;
1867: M2 = tmatd->lda;
1868: for (j=0; j<n; j++) {
1869: for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1870: }
1871: MatDenseRestoreArray(tmat,&v2);
1872: MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);
1873: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1874: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1875: *matout = tmat;
1876: }
1877: return(0);
1878: }
1880: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg)
1881: {
1882: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1883: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1884: PetscInt i;
1885: const PetscScalar *v1,*v2;
1886: PetscErrorCode ierr;
1889: if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1890: if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1891: MatDenseGetArrayRead(A1,&v1);
1892: MatDenseGetArrayRead(A2,&v2);
1893: for (i=0; i<A1->cmap->n; i++) {
1894: PetscArraycmp(v1,v2,A1->rmap->n,flg);
1895: if (*flg == PETSC_FALSE) return(0);
1896: v1 += mat1->lda;
1897: v2 += mat2->lda;
1898: }
1899: MatDenseRestoreArrayRead(A1,&v1);
1900: MatDenseRestoreArrayRead(A2,&v2);
1901: *flg = PETSC_TRUE;
1902: return(0);
1903: }
1905: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1906: {
1907: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1908: PetscInt i,n,len;
1909: PetscScalar *x;
1910: const PetscScalar *vv;
1911: PetscErrorCode ierr;
1914: VecGetSize(v,&n);
1915: VecGetArray(v,&x);
1916: len = PetscMin(A->rmap->n,A->cmap->n);
1917: MatDenseGetArrayRead(A,&vv);
1918: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1919: for (i=0; i<len; i++) {
1920: x[i] = vv[i*mat->lda + i];
1921: }
1922: MatDenseRestoreArrayRead(A,&vv);
1923: VecRestoreArray(v,&x);
1924: return(0);
1925: }
1927: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1928: {
1929: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1930: const PetscScalar *l,*r;
1931: PetscScalar x,*v,*vv;
1932: PetscErrorCode ierr;
1933: PetscInt i,j,m = A->rmap->n,n = A->cmap->n;
1936: MatDenseGetArray(A,&vv);
1937: if (ll) {
1938: VecGetSize(ll,&m);
1939: VecGetArrayRead(ll,&l);
1940: if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1941: for (i=0; i<m; i++) {
1942: x = l[i];
1943: v = vv + i;
1944: for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1945: }
1946: VecRestoreArrayRead(ll,&l);
1947: PetscLogFlops(1.0*n*m);
1948: }
1949: if (rr) {
1950: VecGetSize(rr,&n);
1951: VecGetArrayRead(rr,&r);
1952: if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1953: for (i=0; i<n; i++) {
1954: x = r[i];
1955: v = vv + i*mat->lda;
1956: for (j=0; j<m; j++) (*v++) *= x;
1957: }
1958: VecRestoreArrayRead(rr,&r);
1959: PetscLogFlops(1.0*n*m);
1960: }
1961: MatDenseRestoreArray(A,&vv);
1962: return(0);
1963: }
1965: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1966: {
1967: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1968: PetscScalar *v,*vv;
1969: PetscReal sum = 0.0;
1970: PetscInt lda =mat->lda,m=A->rmap->n,i,j;
1971: PetscErrorCode ierr;
1974: MatDenseGetArrayRead(A,(const PetscScalar**)&vv);
1975: v = vv;
1976: if (type == NORM_FROBENIUS) {
1977: if (lda>m) {
1978: for (j=0; j<A->cmap->n; j++) {
1979: v = vv+j*lda;
1980: for (i=0; i<m; i++) {
1981: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1982: }
1983: }
1984: } else {
1985: #if defined(PETSC_USE_REAL___FP16)
1986: PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1987: PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&cnt,v,&one));
1988: }
1989: #else
1990: for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1991: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1992: }
1993: }
1994: *nrm = PetscSqrtReal(sum);
1995: #endif
1996: PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1997: } else if (type == NORM_1) {
1998: *nrm = 0.0;
1999: for (j=0; j<A->cmap->n; j++) {
2000: v = vv + j*mat->lda;
2001: sum = 0.0;
2002: for (i=0; i<A->rmap->n; i++) {
2003: sum += PetscAbsScalar(*v); v++;
2004: }
2005: if (sum > *nrm) *nrm = sum;
2006: }
2007: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
2008: } else if (type == NORM_INFINITY) {
2009: *nrm = 0.0;
2010: for (j=0; j<A->rmap->n; j++) {
2011: v = vv + j;
2012: sum = 0.0;
2013: for (i=0; i<A->cmap->n; i++) {
2014: sum += PetscAbsScalar(*v); v += mat->lda;
2015: }
2016: if (sum > *nrm) *nrm = sum;
2017: }
2018: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
2019: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
2020: MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);
2021: return(0);
2022: }
2024: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
2025: {
2026: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
2030: switch (op) {
2031: case MAT_ROW_ORIENTED:
2032: aij->roworiented = flg;
2033: break;
2034: case MAT_NEW_NONZERO_LOCATIONS:
2035: case MAT_NEW_NONZERO_LOCATION_ERR:
2036: case MAT_NEW_NONZERO_ALLOCATION_ERR:
2037: case MAT_FORCE_DIAGONAL_ENTRIES:
2038: case MAT_KEEP_NONZERO_PATTERN:
2039: case MAT_IGNORE_OFF_PROC_ENTRIES:
2040: case MAT_USE_HASH_TABLE:
2041: case MAT_IGNORE_ZERO_ENTRIES:
2042: case MAT_IGNORE_LOWER_TRIANGULAR:
2043: case MAT_SORTED_FULL:
2044: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
2045: break;
2046: case MAT_SPD:
2047: case MAT_SYMMETRIC:
2048: case MAT_STRUCTURALLY_SYMMETRIC:
2049: case MAT_HERMITIAN:
2050: case MAT_SYMMETRY_ETERNAL:
2051: /* These options are handled directly by MatSetOption() */
2052: break;
2053: default:
2054: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
2055: }
2056: return(0);
2057: }
2059: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
2060: {
2061: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
2063: PetscInt lda=l->lda,m=A->rmap->n,n=A->cmap->n,j;
2064: PetscScalar *v;
2067: MatDenseGetArrayWrite(A,&v);
2068: if (lda>m) {
2069: for (j=0; j<n; j++) {
2070: PetscArrayzero(v+j*lda,m);
2071: }
2072: } else {
2073: PetscArrayzero(v,PetscInt64Mult(m,n));
2074: }
2075: MatDenseRestoreArrayWrite(A,&v);
2076: return(0);
2077: }
2079: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
2080: {
2081: PetscErrorCode ierr;
2082: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
2083: PetscInt m = l->lda, n = A->cmap->n, i,j;
2084: PetscScalar *slot,*bb,*v;
2085: const PetscScalar *xx;
2088: if (PetscDefined(USE_DEBUG)) {
2089: for (i=0; i<N; i++) {
2090: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
2091: 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);
2092: }
2093: }
2094: if (!N) return(0);
2096: /* fix right hand side if needed */
2097: if (x && b) {
2098: VecGetArrayRead(x,&xx);
2099: VecGetArray(b,&bb);
2100: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
2101: VecRestoreArrayRead(x,&xx);
2102: VecRestoreArray(b,&bb);
2103: }
2105: MatDenseGetArray(A,&v);
2106: for (i=0; i<N; i++) {
2107: slot = v + rows[i];
2108: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
2109: }
2110: if (diag != 0.0) {
2111: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
2112: for (i=0; i<N; i++) {
2113: slot = v + (m+1)*rows[i];
2114: *slot = diag;
2115: }
2116: }
2117: MatDenseRestoreArray(A,&v);
2118: return(0);
2119: }
2121: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda)
2122: {
2123: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
2126: *lda = mat->lda;
2127: return(0);
2128: }
2130: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array)
2131: {
2132: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
2135: if (mat->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
2136: *array = mat->v;
2137: return(0);
2138: }
2140: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array)
2141: {
2143: *array = NULL;
2144: return(0);
2145: }
2147: /*@
2148: MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray()
2150: Not collective
2152: Input Parameter:
2153: . mat - a MATSEQDENSE or MATMPIDENSE matrix
2155: Output Parameter:
2156: . lda - the leading dimension
2158: Level: intermediate
2160: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseSetLDA()
2161: @*/
2162: PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda)
2163: {
2169: PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));
2170: return(0);
2171: }
2173: /*@
2174: MatDenseSetLDA - Sets the leading dimension of the array used by the dense matrix
2176: Not collective
2178: Input Parameters:
2179: + mat - a MATSEQDENSE or MATMPIDENSE matrix
2180: - lda - the leading dimension
2182: Level: intermediate
2184: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetLDA()
2185: @*/
2186: PetscErrorCode MatDenseSetLDA(Mat A,PetscInt lda)
2187: {
2192: PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda));
2193: return(0);
2194: }
2196: /*@C
2197: MatDenseGetArray - gives read-write access to the array where the data for a dense matrix is stored
2199: Logically Collective on Mat
2201: Input Parameter:
2202: . mat - a dense matrix
2204: Output Parameter:
2205: . array - pointer to the data
2207: Level: intermediate
2209: .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2210: @*/
2211: PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array)
2212: {
2218: PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
2219: return(0);
2220: }
2222: /*@C
2223: MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
2225: Logically Collective on Mat
2227: Input Parameters:
2228: + mat - a dense matrix
2229: - array - pointer to the data
2231: Level: intermediate
2233: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2234: @*/
2235: PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array)
2236: {
2242: PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
2243: PetscObjectStateIncrease((PetscObject)A);
2244: #if defined(PETSC_HAVE_CUDA)
2245: A->offloadmask = PETSC_OFFLOAD_CPU;
2246: #endif
2247: return(0);
2248: }
2250: /*@C
2251: MatDenseGetArrayRead - gives read-only access to the array where the data for a dense matrix is stored
2253: Not Collective
2255: Input Parameter:
2256: . mat - a dense matrix
2258: Output Parameter:
2259: . array - pointer to the data
2261: Level: intermediate
2263: .seealso: MatDenseRestoreArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2264: @*/
2265: PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array)
2266: {
2272: PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
2273: return(0);
2274: }
2276: /*@C
2277: MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayRead()
2279: Not Collective
2281: Input Parameters:
2282: + mat - a dense matrix
2283: - array - pointer to the data
2285: Level: intermediate
2287: .seealso: MatDenseGetArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
2288: @*/
2289: PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
2290: {
2296: PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
2297: return(0);
2298: }
2300: /*@C
2301: MatDenseGetArrayWrite - gives write-only access to the array where the data for a dense matrix is stored
2303: Not Collective
2305: Input Parameter:
2306: . mat - a dense matrix
2308: Output Parameter:
2309: . array - pointer to the data
2311: Level: intermediate
2313: .seealso: MatDenseRestoreArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
2314: @*/
2315: PetscErrorCode MatDenseGetArrayWrite(Mat A,PetscScalar **array)
2316: {
2322: PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array));
2323: return(0);
2324: }
2326: /*@C
2327: MatDenseRestoreArrayWrite - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayWrite()
2329: Not Collective
2331: Input Parameters:
2332: + mat - a dense matrix
2333: - array - pointer to the data
2335: Level: intermediate
2337: .seealso: MatDenseGetArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
2338: @*/
2339: PetscErrorCode MatDenseRestoreArrayWrite(Mat A,PetscScalar **array)
2340: {
2346: PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array));
2347: PetscObjectStateIncrease((PetscObject)A);
2348: #if defined(PETSC_HAVE_CUDA)
2349: A->offloadmask = PETSC_OFFLOAD_CPU;
2350: #endif
2351: return(0);
2352: }
2354: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B)
2355: {
2356: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
2358: PetscInt i,j,nrows,ncols,ldb;
2359: const PetscInt *irow,*icol;
2360: PetscScalar *av,*bv,*v = mat->v;
2361: Mat newmat;
2364: ISGetIndices(isrow,&irow);
2365: ISGetIndices(iscol,&icol);
2366: ISGetLocalSize(isrow,&nrows);
2367: ISGetLocalSize(iscol,&ncols);
2369: /* Check submatrixcall */
2370: if (scall == MAT_REUSE_MATRIX) {
2371: PetscInt n_cols,n_rows;
2372: MatGetSize(*B,&n_rows,&n_cols);
2373: if (n_rows != nrows || n_cols != ncols) {
2374: /* resize the result matrix to match number of requested rows/columns */
2375: MatSetSizes(*B,nrows,ncols,nrows,ncols);
2376: }
2377: newmat = *B;
2378: } else {
2379: /* Create and fill new matrix */
2380: MatCreate(PetscObjectComm((PetscObject)A),&newmat);
2381: MatSetSizes(newmat,nrows,ncols,nrows,ncols);
2382: MatSetType(newmat,((PetscObject)A)->type_name);
2383: MatSeqDenseSetPreallocation(newmat,NULL);
2384: }
2386: /* Now extract the data pointers and do the copy,column at a time */
2387: MatDenseGetArray(newmat,&bv);
2388: MatDenseGetLDA(newmat,&ldb);
2389: for (i=0; i<ncols; i++) {
2390: av = v + mat->lda*icol[i];
2391: for (j=0; j<nrows; j++) bv[j] = av[irow[j]];
2392: bv += ldb;
2393: }
2394: MatDenseRestoreArray(newmat,&bv);
2396: /* Assemble the matrices so that the correct flags are set */
2397: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
2398: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
2400: /* Free work space */
2401: ISRestoreIndices(isrow,&irow);
2402: ISRestoreIndices(iscol,&icol);
2403: *B = newmat;
2404: return(0);
2405: }
2407: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
2408: {
2410: PetscInt i;
2413: if (scall == MAT_INITIAL_MATRIX) {
2414: PetscCalloc1(n,B);
2415: }
2417: for (i=0; i<n; i++) {
2418: MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],scall,&(*B)[i]);
2419: }
2420: return(0);
2421: }
2423: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
2424: {
2426: return(0);
2427: }
2429: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
2430: {
2432: return(0);
2433: }
2435: PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
2436: {
2437: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
2438: PetscErrorCode ierr;
2439: const PetscScalar *va;
2440: PetscScalar *vb;
2441: PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;
2444: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2445: if (A->ops->copy != B->ops->copy) {
2446: MatCopy_Basic(A,B,str);
2447: return(0);
2448: }
2449: if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
2450: MatDenseGetArrayRead(A,&va);
2451: MatDenseGetArray(B,&vb);
2452: if (lda1>m || lda2>m) {
2453: for (j=0; j<n; j++) {
2454: PetscArraycpy(vb+j*lda2,va+j*lda1,m);
2455: }
2456: } else {
2457: PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);
2458: }
2459: MatDenseRestoreArray(B,&vb);
2460: MatDenseRestoreArrayRead(A,&va);
2461: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2462: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2463: return(0);
2464: }
2466: static PetscErrorCode MatSetUp_SeqDense(Mat A)
2467: {
2471: PetscLayoutSetUp(A->rmap);
2472: PetscLayoutSetUp(A->cmap);
2473: if (!A->preallocated) {
2474: MatSeqDenseSetPreallocation(A,NULL);
2475: }
2476: return(0);
2477: }
2479: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2480: {
2481: Mat_SeqDense *mat = (Mat_SeqDense *) A->data;
2482: PetscInt i,nz = A->rmap->n*A->cmap->n;
2483: PetscInt min = PetscMin(A->rmap->n,A->cmap->n);
2484: PetscScalar *aa;
2488: MatDenseGetArray(A,&aa);
2489: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
2490: MatDenseRestoreArray(A,&aa);
2491: if (mat->tau) for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]);
2492: return(0);
2493: }
2495: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2496: {
2497: PetscInt i,nz = A->rmap->n*A->cmap->n;
2498: PetscScalar *aa;
2502: MatDenseGetArray(A,&aa);
2503: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
2504: MatDenseRestoreArray(A,&aa);
2505: return(0);
2506: }
2508: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2509: {
2510: PetscInt i,nz = A->rmap->n*A->cmap->n;
2511: PetscScalar *aa;
2515: MatDenseGetArray(A,&aa);
2516: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2517: MatDenseRestoreArray(A,&aa);
2518: return(0);
2519: }
2521: /* ----------------------------------------------------------------*/
2522: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2523: {
2525: PetscInt m=A->rmap->n,n=B->cmap->n;
2526: PetscBool cisdense;
2529: MatSetSizes(C,m,n,m,n);
2530: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2531: if (!cisdense) {
2532: PetscBool flg;
2534: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2535: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2536: }
2537: MatSetUp(C);
2538: return(0);
2539: }
2541: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2542: {
2543: Mat_SeqDense *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data;
2544: PetscBLASInt m,n,k;
2545: const PetscScalar *av,*bv;
2546: PetscScalar *cv;
2547: PetscScalar _DOne=1.0,_DZero=0.0;
2548: PetscErrorCode ierr;
2551: PetscBLASIntCast(C->rmap->n,&m);
2552: PetscBLASIntCast(C->cmap->n,&n);
2553: PetscBLASIntCast(A->cmap->n,&k);
2554: if (!m || !n || !k) return(0);
2555: MatDenseGetArrayRead(A,&av);
2556: MatDenseGetArrayRead(B,&bv);
2557: MatDenseGetArrayWrite(C,&cv);
2558: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2559: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2560: MatDenseRestoreArrayRead(A,&av);
2561: MatDenseRestoreArrayRead(B,&bv);
2562: MatDenseRestoreArrayWrite(C,&cv);
2563: return(0);
2564: }
2566: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2567: {
2569: PetscInt m=A->rmap->n,n=B->rmap->n;
2570: PetscBool cisdense;
2573: MatSetSizes(C,m,n,m,n);
2574: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2575: if (!cisdense) {
2576: PetscBool flg;
2578: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2579: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2580: }
2581: MatSetUp(C);
2582: return(0);
2583: }
2585: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2586: {
2587: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2588: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2589: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2590: const PetscScalar *av,*bv;
2591: PetscScalar *cv;
2592: PetscBLASInt m,n,k;
2593: PetscScalar _DOne=1.0,_DZero=0.0;
2594: PetscErrorCode ierr;
2597: PetscBLASIntCast(C->rmap->n,&m);
2598: PetscBLASIntCast(C->cmap->n,&n);
2599: PetscBLASIntCast(A->cmap->n,&k);
2600: if (!m || !n || !k) return(0);
2601: MatDenseGetArrayRead(A,&av);
2602: MatDenseGetArrayRead(B,&bv);
2603: MatDenseGetArrayWrite(C,&cv);
2604: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2605: MatDenseRestoreArrayRead(A,&av);
2606: MatDenseRestoreArrayRead(B,&bv);
2607: MatDenseRestoreArrayWrite(C,&cv);
2608: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2609: return(0);
2610: }
2612: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2613: {
2615: PetscInt m=A->cmap->n,n=B->cmap->n;
2616: PetscBool cisdense;
2619: MatSetSizes(C,m,n,m,n);
2620: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2621: if (!cisdense) {
2622: PetscBool flg;
2624: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2625: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2626: }
2627: MatSetUp(C);
2628: return(0);
2629: }
2631: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2632: {
2633: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2634: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2635: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2636: const PetscScalar *av,*bv;
2637: PetscScalar *cv;
2638: PetscBLASInt m,n,k;
2639: PetscScalar _DOne=1.0,_DZero=0.0;
2640: PetscErrorCode ierr;
2643: PetscBLASIntCast(C->rmap->n,&m);
2644: PetscBLASIntCast(C->cmap->n,&n);
2645: PetscBLASIntCast(A->rmap->n,&k);
2646: if (!m || !n || !k) return(0);
2647: MatDenseGetArrayRead(A,&av);
2648: MatDenseGetArrayRead(B,&bv);
2649: MatDenseGetArrayWrite(C,&cv);
2650: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2651: MatDenseRestoreArrayRead(A,&av);
2652: MatDenseRestoreArrayRead(B,&bv);
2653: MatDenseRestoreArrayWrite(C,&cv);
2654: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2655: return(0);
2656: }
2658: /* ----------------------------------------------- */
2659: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2660: {
2662: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2663: C->ops->productsymbolic = MatProductSymbolic_AB;
2664: return(0);
2665: }
2667: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2668: {
2670: C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2671: C->ops->productsymbolic = MatProductSymbolic_AtB;
2672: return(0);
2673: }
2675: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2676: {
2678: C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2679: C->ops->productsymbolic = MatProductSymbolic_ABt;
2680: return(0);
2681: }
2683: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2684: {
2686: Mat_Product *product = C->product;
2689: switch (product->type) {
2690: case MATPRODUCT_AB:
2691: MatProductSetFromOptions_SeqDense_AB(C);
2692: break;
2693: case MATPRODUCT_AtB:
2694: MatProductSetFromOptions_SeqDense_AtB(C);
2695: break;
2696: case MATPRODUCT_ABt:
2697: MatProductSetFromOptions_SeqDense_ABt(C);
2698: break;
2699: default:
2700: break;
2701: }
2702: return(0);
2703: }
2704: /* ----------------------------------------------- */
2706: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2707: {
2708: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2709: PetscErrorCode ierr;
2710: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2711: PetscScalar *x;
2712: const PetscScalar *aa;
2715: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2716: VecGetArray(v,&x);
2717: VecGetLocalSize(v,&p);
2718: MatDenseGetArrayRead(A,&aa);
2719: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2720: for (i=0; i<m; i++) {
2721: x[i] = aa[i]; if (idx) idx[i] = 0;
2722: for (j=1; j<n; j++) {
2723: if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2724: }
2725: }
2726: MatDenseRestoreArrayRead(A,&aa);
2727: VecRestoreArray(v,&x);
2728: return(0);
2729: }
2731: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2732: {
2733: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2734: PetscErrorCode ierr;
2735: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2736: PetscScalar *x;
2737: PetscReal atmp;
2738: const PetscScalar *aa;
2741: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2742: VecGetArray(v,&x);
2743: VecGetLocalSize(v,&p);
2744: MatDenseGetArrayRead(A,&aa);
2745: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2746: for (i=0; i<m; i++) {
2747: x[i] = PetscAbsScalar(aa[i]);
2748: for (j=1; j<n; j++) {
2749: atmp = PetscAbsScalar(aa[i+a->lda*j]);
2750: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2751: }
2752: }
2753: MatDenseRestoreArrayRead(A,&aa);
2754: VecRestoreArray(v,&x);
2755: return(0);
2756: }
2758: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2759: {
2760: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2761: PetscErrorCode ierr;
2762: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2763: PetscScalar *x;
2764: const PetscScalar *aa;
2767: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2768: MatDenseGetArrayRead(A,&aa);
2769: VecGetArray(v,&x);
2770: VecGetLocalSize(v,&p);
2771: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2772: for (i=0; i<m; i++) {
2773: x[i] = aa[i]; if (idx) idx[i] = 0;
2774: for (j=1; j<n; j++) {
2775: if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2776: }
2777: }
2778: VecRestoreArray(v,&x);
2779: MatDenseRestoreArrayRead(A,&aa);
2780: return(0);
2781: }
2783: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2784: {
2785: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2786: PetscErrorCode ierr;
2787: PetscScalar *x;
2788: const PetscScalar *aa;
2791: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2792: MatDenseGetArrayRead(A,&aa);
2793: VecGetArray(v,&x);
2794: PetscArraycpy(x,aa+col*a->lda,A->rmap->n);
2795: VecRestoreArray(v,&x);
2796: MatDenseRestoreArrayRead(A,&aa);
2797: return(0);
2798: }
2800: PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A,PetscInt type,PetscReal *reductions)
2801: {
2802: PetscErrorCode ierr;
2803: PetscInt i,j,m,n;
2804: const PetscScalar *a;
2807: MatGetSize(A,&m,&n);
2808: PetscArrayzero(reductions,n);
2809: MatDenseGetArrayRead(A,&a);
2810: if (type == NORM_2) {
2811: for (i=0; i<n; i++) {
2812: for (j=0; j<m; j++) {
2813: reductions[i] += PetscAbsScalar(a[j]*a[j]);
2814: }
2815: a += m;
2816: }
2817: } else if (type == NORM_1) {
2818: for (i=0; i<n; i++) {
2819: for (j=0; j<m; j++) {
2820: reductions[i] += PetscAbsScalar(a[j]);
2821: }
2822: a += m;
2823: }
2824: } else if (type == NORM_INFINITY) {
2825: for (i=0; i<n; i++) {
2826: for (j=0; j<m; j++) {
2827: reductions[i] = PetscMax(PetscAbsScalar(a[j]),reductions[i]);
2828: }
2829: a += m;
2830: }
2831: } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
2832: for (i=0; i<n; i++) {
2833: for (j=0; j<m; j++) {
2834: reductions[i] += PetscRealPart(a[j]);
2835: }
2836: a += m;
2837: }
2838: } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
2839: for (i=0; i<n; i++) {
2840: for (j=0; j<m; j++) {
2841: reductions[i] += PetscImaginaryPart(a[j]);
2842: }
2843: a += m;
2844: }
2845: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown reduction type");
2846: MatDenseRestoreArrayRead(A,&a);
2847: if (type == NORM_2) {
2848: for (i=0; i<n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
2849: } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
2850: for (i=0; i<n; i++) reductions[i] /= m;
2851: }
2852: return(0);
2853: }
2855: static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2856: {
2858: PetscScalar *a;
2859: PetscInt lda,m,n,i,j;
2862: MatGetSize(x,&m,&n);
2863: MatDenseGetLDA(x,&lda);
2864: MatDenseGetArray(x,&a);
2865: for (j=0; j<n; j++) {
2866: for (i=0; i<m; i++) {
2867: PetscRandomGetValue(rctx,a+j*lda+i);
2868: }
2869: }
2870: MatDenseRestoreArray(x,&a);
2871: return(0);
2872: }
2874: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d)
2875: {
2877: *missing = PETSC_FALSE;
2878: return(0);
2879: }
2881: /* vals is not const */
2882: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2883: {
2885: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2886: PetscScalar *v;
2889: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2890: MatDenseGetArray(A,&v);
2891: *vals = v+col*a->lda;
2892: MatDenseRestoreArray(A,&v);
2893: return(0);
2894: }
2896: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2897: {
2899: *vals = NULL; /* user cannot accidentally use the array later */
2900: return(0);
2901: }
2903: /* -------------------------------------------------------------------*/
2904: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2905: MatGetRow_SeqDense,
2906: MatRestoreRow_SeqDense,
2907: MatMult_SeqDense,
2908: /* 4*/ MatMultAdd_SeqDense,
2909: MatMultTranspose_SeqDense,
2910: MatMultTransposeAdd_SeqDense,
2911: NULL,
2912: NULL,
2913: NULL,
2914: /* 10*/ NULL,
2915: MatLUFactor_SeqDense,
2916: MatCholeskyFactor_SeqDense,
2917: MatSOR_SeqDense,
2918: MatTranspose_SeqDense,
2919: /* 15*/ MatGetInfo_SeqDense,
2920: MatEqual_SeqDense,
2921: MatGetDiagonal_SeqDense,
2922: MatDiagonalScale_SeqDense,
2923: MatNorm_SeqDense,
2924: /* 20*/ MatAssemblyBegin_SeqDense,
2925: MatAssemblyEnd_SeqDense,
2926: MatSetOption_SeqDense,
2927: MatZeroEntries_SeqDense,
2928: /* 24*/ MatZeroRows_SeqDense,
2929: NULL,
2930: NULL,
2931: NULL,
2932: NULL,
2933: /* 29*/ MatSetUp_SeqDense,
2934: NULL,
2935: NULL,
2936: NULL,
2937: NULL,
2938: /* 34*/ MatDuplicate_SeqDense,
2939: NULL,
2940: NULL,
2941: NULL,
2942: NULL,
2943: /* 39*/ MatAXPY_SeqDense,
2944: MatCreateSubMatrices_SeqDense,
2945: NULL,
2946: MatGetValues_SeqDense,
2947: MatCopy_SeqDense,
2948: /* 44*/ MatGetRowMax_SeqDense,
2949: MatScale_SeqDense,
2950: MatShift_Basic,
2951: NULL,
2952: MatZeroRowsColumns_SeqDense,
2953: /* 49*/ MatSetRandom_SeqDense,
2954: NULL,
2955: NULL,
2956: NULL,
2957: NULL,
2958: /* 54*/ NULL,
2959: NULL,
2960: NULL,
2961: NULL,
2962: NULL,
2963: /* 59*/ MatCreateSubMatrix_SeqDense,
2964: MatDestroy_SeqDense,
2965: MatView_SeqDense,
2966: NULL,
2967: NULL,
2968: /* 64*/ NULL,
2969: NULL,
2970: NULL,
2971: NULL,
2972: NULL,
2973: /* 69*/ MatGetRowMaxAbs_SeqDense,
2974: NULL,
2975: NULL,
2976: NULL,
2977: NULL,
2978: /* 74*/ NULL,
2979: NULL,
2980: NULL,
2981: NULL,
2982: NULL,
2983: /* 79*/ NULL,
2984: NULL,
2985: NULL,
2986: NULL,
2987: /* 83*/ MatLoad_SeqDense,
2988: MatIsSymmetric_SeqDense,
2989: MatIsHermitian_SeqDense,
2990: NULL,
2991: NULL,
2992: NULL,
2993: /* 89*/ NULL,
2994: NULL,
2995: MatMatMultNumeric_SeqDense_SeqDense,
2996: NULL,
2997: NULL,
2998: /* 94*/ NULL,
2999: NULL,
3000: NULL,
3001: MatMatTransposeMultNumeric_SeqDense_SeqDense,
3002: NULL,
3003: /* 99*/ MatProductSetFromOptions_SeqDense,
3004: NULL,
3005: NULL,
3006: MatConjugate_SeqDense,
3007: NULL,
3008: /*104*/ NULL,
3009: MatRealPart_SeqDense,
3010: MatImaginaryPart_SeqDense,
3011: NULL,
3012: NULL,
3013: /*109*/ NULL,
3014: NULL,
3015: MatGetRowMin_SeqDense,
3016: MatGetColumnVector_SeqDense,
3017: MatMissingDiagonal_SeqDense,
3018: /*114*/ NULL,
3019: NULL,
3020: NULL,
3021: NULL,
3022: NULL,
3023: /*119*/ NULL,
3024: NULL,
3025: NULL,
3026: NULL,
3027: NULL,
3028: /*124*/ NULL,
3029: MatGetColumnReductions_SeqDense,
3030: NULL,
3031: NULL,
3032: NULL,
3033: /*129*/ NULL,
3034: NULL,
3035: NULL,
3036: MatTransposeMatMultNumeric_SeqDense_SeqDense,
3037: NULL,
3038: /*134*/ NULL,
3039: NULL,
3040: NULL,
3041: NULL,
3042: NULL,
3043: /*139*/ NULL,
3044: NULL,
3045: NULL,
3046: NULL,
3047: NULL,
3048: MatCreateMPIMatConcatenateSeqMat_SeqDense,
3049: /*145*/ NULL,
3050: NULL,
3051: NULL
3052: };
3054: /*@C
3055: MatCreateSeqDense - Creates a sequential dense matrix that
3056: is stored in column major order (the usual Fortran 77 manner). Many
3057: of the matrix operations use the BLAS and LAPACK routines.
3059: Collective
3061: Input Parameters:
3062: + comm - MPI communicator, set to PETSC_COMM_SELF
3063: . m - number of rows
3064: . n - number of columns
3065: - data - optional location of matrix data in column major order. Set data=NULL for PETSc
3066: to control all matrix memory allocation.
3068: Output Parameter:
3069: . A - the matrix
3071: Notes:
3072: The data input variable is intended primarily for Fortran programmers
3073: who wish to allocate their own matrix memory space. Most users should
3074: set data=NULL.
3076: Level: intermediate
3078: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
3079: @*/
3080: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
3081: {
3085: MatCreate(comm,A);
3086: MatSetSizes(*A,m,n,m,n);
3087: MatSetType(*A,MATSEQDENSE);
3088: MatSeqDenseSetPreallocation(*A,data);
3089: return(0);
3090: }
3092: /*@C
3093: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
3095: Collective
3097: Input Parameters:
3098: + B - the matrix
3099: - data - the array (or NULL)
3101: Notes:
3102: The data input variable is intended primarily for Fortran programmers
3103: who wish to allocate their own matrix memory space. Most users should
3104: need not call this routine.
3106: Level: intermediate
3108: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatDenseSetLDA()
3110: @*/
3111: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
3112: {
3117: PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
3118: return(0);
3119: }
3121: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
3122: {
3123: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
3127: if (b->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3128: B->preallocated = PETSC_TRUE;
3130: PetscLayoutSetUp(B->rmap);
3131: PetscLayoutSetUp(B->cmap);
3133: if (b->lda <= 0) b->lda = B->rmap->n;
3135: if (!data) { /* petsc-allocated storage */
3136: if (!b->user_alloc) { PetscFree(b->v); }
3137: PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v);
3138: PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar));
3140: b->user_alloc = PETSC_FALSE;
3141: } else { /* user-allocated storage */
3142: if (!b->user_alloc) { PetscFree(b->v); }
3143: b->v = data;
3144: b->user_alloc = PETSC_TRUE;
3145: }
3146: B->assembled = PETSC_TRUE;
3147: return(0);
3148: }
3150: #if defined(PETSC_HAVE_ELEMENTAL)
3151: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
3152: {
3153: Mat mat_elemental;
3154: PetscErrorCode ierr;
3155: const PetscScalar *array;
3156: PetscScalar *v_colwise;
3157: PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;
3160: PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
3161: MatDenseGetArrayRead(A,&array);
3162: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
3163: k = 0;
3164: for (j=0; j<N; j++) {
3165: cols[j] = j;
3166: for (i=0; i<M; i++) {
3167: v_colwise[j*M+i] = array[k++];
3168: }
3169: }
3170: for (i=0; i<M; i++) {
3171: rows[i] = i;
3172: }
3173: MatDenseRestoreArrayRead(A,&array);
3175: MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
3176: MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
3177: MatSetType(mat_elemental,MATELEMENTAL);
3178: MatSetUp(mat_elemental);
3180: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
3181: MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
3182: MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
3183: MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
3184: PetscFree3(v_colwise,rows,cols);
3186: if (reuse == MAT_INPLACE_MATRIX) {
3187: MatHeaderReplace(A,&mat_elemental);
3188: } else {
3189: *newmat = mat_elemental;
3190: }
3191: return(0);
3192: }
3193: #endif
3195: PetscErrorCode MatDenseSetLDA_SeqDense(Mat B,PetscInt lda)
3196: {
3197: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
3198: PetscBool data;
3201: data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
3202: if (!b->user_alloc && data && b->lda!=lda) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"LDA cannot be changed after allocation of internal storage");
3203: 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);
3204: b->lda = lda;
3205: return(0);
3206: }
3208: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
3209: {
3211: PetscMPIInt size;
3214: MPI_Comm_size(comm,&size);
3215: if (size == 1) {
3216: if (scall == MAT_INITIAL_MATRIX) {
3217: MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
3218: } else {
3219: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
3220: }
3221: } else {
3222: MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
3223: }
3224: return(0);
3225: }
3227: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
3228: {
3229: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3233: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
3234: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3235: if (!a->cvec) {
3236: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3237: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3238: }
3239: a->vecinuse = col + 1;
3240: MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);
3241: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3242: *v = a->cvec;
3243: return(0);
3244: }
3246: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
3247: {
3248: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3252: if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
3253: if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
3254: a->vecinuse = 0;
3255: MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);
3256: VecResetArray(a->cvec);
3257: *v = NULL;
3258: return(0);
3259: }
3261: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
3262: {
3263: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3267: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
3268: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3269: if (!a->cvec) {
3270: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3271: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3272: }
3273: a->vecinuse = col + 1;
3274: MatDenseGetArrayRead(A,&a->ptrinuse);
3275: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3276: VecLockReadPush(a->cvec);
3277: *v = a->cvec;
3278: return(0);
3279: }
3281: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
3282: {
3283: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3287: if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
3288: if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
3289: a->vecinuse = 0;
3290: MatDenseRestoreArrayRead(A,&a->ptrinuse);
3291: VecLockReadPop(a->cvec);
3292: VecResetArray(a->cvec);
3293: *v = NULL;
3294: return(0);
3295: }
3297: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
3298: {
3299: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3303: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
3304: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3305: if (!a->cvec) {
3306: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
3307: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
3308: }
3309: a->vecinuse = col + 1;
3310: MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);
3311: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
3312: *v = a->cvec;
3313: return(0);
3314: }
3316: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
3317: {
3318: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3322: if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
3323: if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
3324: a->vecinuse = 0;
3325: MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);
3326: VecResetArray(a->cvec);
3327: *v = NULL;
3328: return(0);
3329: }
3331: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3332: {
3333: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3337: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
3338: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
3339: if (a->cmat && cend-cbegin != a->cmat->cmap->N) {
3340: MatDestroy(&a->cmat);
3341: }
3342: if (!a->cmat) {
3343: 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);
3344: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat);
3345: } else {
3346: MatDensePlaceArray(a->cmat,a->v + (size_t)cbegin * (size_t)a->lda);
3347: }
3348: MatDenseSetLDA(a->cmat,a->lda);
3349: a->matinuse = cbegin + 1;
3350: *v = a->cmat;
3351: return(0);
3352: }
3354: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v)
3355: {
3356: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3360: if (!a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetSubMatrix() first");
3361: if (!a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column matrix");
3362: if (*v != a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not the matrix obtained from MatDenseGetSubMatrix()");
3363: a->matinuse = 0;
3364: MatDenseResetArray(a->cmat);
3365: *v = NULL;
3366: return(0);
3367: }
3369: /*MC
3370: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
3372: Options Database Keys:
3373: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
3375: Level: beginner
3377: .seealso: MatCreateSeqDense()
3379: M*/
3380: PetscErrorCode MatCreate_SeqDense(Mat B)
3381: {
3382: Mat_SeqDense *b;
3384: PetscMPIInt size;
3387: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
3388: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
3390: PetscNewLog(B,&b);
3391: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
3392: B->data = (void*)b;
3394: b->roworiented = PETSC_TRUE;
3396: PetscObjectComposeFunction((PetscObject)B,"MatQRFactor_C",MatQRFactor_SeqDense);
3397: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);
3398: PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense);
3399: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
3400: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
3401: PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
3402: PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
3403: PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense);
3404: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
3405: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);
3406: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense);
3407: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense);
3408: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
3409: #if defined(PETSC_HAVE_ELEMENTAL)
3410: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
3411: #endif
3412: #if defined(PETSC_HAVE_SCALAPACK)
3413: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_scalapack_C",MatConvert_Dense_ScaLAPACK);
3414: #endif
3415: #if defined(PETSC_HAVE_CUDA)
3416: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);
3417: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3418: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);
3419: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3420: #endif
3421: PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
3422: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);
3423: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);
3424: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
3425: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
3427: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
3428: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
3429: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);
3430: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);
3431: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);
3432: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);
3433: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);
3434: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);
3435: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense);
3436: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense);
3437: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
3438: return(0);
3439: }
3441: /*@C
3442: 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.
3444: Not Collective
3446: Input Parameters:
3447: + mat - a MATSEQDENSE or MATMPIDENSE matrix
3448: - col - column index
3450: Output Parameter:
3451: . vals - pointer to the data
3453: Level: intermediate
3455: .seealso: MatDenseRestoreColumn()
3456: @*/
3457: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
3458: {
3465: PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
3466: return(0);
3467: }
3469: /*@C
3470: MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn().
3472: Not Collective
3474: Input Parameter:
3475: . mat - a MATSEQDENSE or MATMPIDENSE matrix
3477: Output Parameter:
3478: . vals - pointer to the data
3480: Level: intermediate
3482: .seealso: MatDenseGetColumn()
3483: @*/
3484: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
3485: {
3491: PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
3492: return(0);
3493: }
3495: /*@
3496: MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec.
3498: Collective
3500: Input Parameters:
3501: + mat - the Mat object
3502: - col - the column index
3504: Output Parameter:
3505: . v - the vector
3507: Notes:
3508: The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVec() when the vector is no longer needed.
3509: Use MatDenseGetColumnVecRead() to obtain read-only access or MatDenseGetColumnVecWrite() for write-only access.
3511: Level: intermediate
3513: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3514: @*/
3515: PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v)
3516: {
3524: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3525: 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);
3526: PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3527: return(0);
3528: }
3530: /*@
3531: MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec().
3533: Collective
3535: Input Parameters:
3536: + mat - the Mat object
3537: . col - the column index
3538: - v - the Vec object
3540: Level: intermediate
3542: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3543: @*/
3544: PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v)
3545: {
3553: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3554: 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);
3555: PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3556: return(0);
3557: }
3559: /*@
3560: MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec.
3562: Collective
3564: Input Parameters:
3565: + mat - the Mat object
3566: - col - the column index
3568: Output Parameter:
3569: . v - the vector
3571: Notes:
3572: The vector is owned by PETSc and users cannot modify it.
3573: Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed.
3574: Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access.
3576: Level: intermediate
3578: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3579: @*/
3580: PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v)
3581: {
3589: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3590: 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);
3591: PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3592: return(0);
3593: }
3595: /*@
3596: MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead().
3598: Collective
3600: Input Parameters:
3601: + mat - the Mat object
3602: . col - the column index
3603: - v - the Vec object
3605: Level: intermediate
3607: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite()
3608: @*/
3609: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v)
3610: {
3618: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3619: 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);
3620: PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3621: return(0);
3622: }
3624: /*@
3625: MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec.
3627: Collective
3629: Input Parameters:
3630: + mat - the Mat object
3631: - col - the column index
3633: Output Parameter:
3634: . v - the vector
3636: Notes:
3637: The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed.
3638: Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access.
3640: Level: intermediate
3642: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3643: @*/
3644: PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v)
3645: {
3653: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3654: 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);
3655: PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3656: return(0);
3657: }
3659: /*@
3660: MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite().
3662: Collective
3664: Input Parameters:
3665: + mat - the Mat object
3666: . col - the column index
3667: - v - the Vec object
3669: Level: intermediate
3671: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead()
3672: @*/
3673: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v)
3674: {
3682: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3683: 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);
3684: PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3685: return(0);
3686: }
3688: /*@
3689: MatDenseGetSubMatrix - Gives access to a block of columns of a dense matrix, represented as a Mat.
3691: Collective
3693: Input Parameters:
3694: + mat - the Mat object
3695: . cbegin - the first index in the block
3696: - cend - the last index in the block
3698: Output Parameter:
3699: . v - the matrix
3701: Notes:
3702: The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed.
3704: Level: intermediate
3706: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseRestoreSubMatrix()
3707: @*/
3708: PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3709: {
3718: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3719: 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);
3720: 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);
3721: PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));
3722: return(0);
3723: }
3725: /*@
3726: MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix().
3728: Collective
3730: Input Parameters:
3731: + mat - the Mat object
3732: - v - the Mat object
3734: Level: intermediate
3736: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseGetSubMatrix()
3737: @*/
3738: PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v)
3739: {
3746: PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v));
3747: return(0);
3748: }