Actual source code: dense.c
petsc-3.14.6 2021-03-30
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;
194: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
195: MatScalar *av=a->a;
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: for (i=0; i<m; i++) {
214: PetscInt j;
215: for (j=0;j<ai[1]-ai[0];j++) {
216: b->v[*aj*m+i] = *av;
217: aj++;
218: av++;
219: }
220: ai++;
221: }
223: if (reuse == MAT_INPLACE_MATRIX) {
224: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
225: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
226: MatHeaderReplace(A,&B);
227: } else {
228: if (B) *newmat = B;
229: MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
230: MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
231: }
232: return(0);
233: }
235: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
236: {
237: Mat B;
238: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
240: PetscInt i, j;
241: PetscInt *rows, *nnz;
242: MatScalar *aa = a->v, *vals;
245: MatCreate(PetscObjectComm((PetscObject)A),&B);
246: MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
247: MatSetType(B,MATSEQAIJ);
248: PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);
249: for (j=0; j<A->cmap->n; j++) {
250: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i];
251: aa += a->lda;
252: }
253: MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);
254: aa = a->v;
255: for (j=0; j<A->cmap->n; j++) {
256: PetscInt numRows = 0;
257: for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];}
258: MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);
259: aa += a->lda;
260: }
261: PetscFree3(rows,nnz,vals);
262: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
263: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
265: if (reuse == MAT_INPLACE_MATRIX) {
266: MatHeaderReplace(A,&B);
267: } else {
268: *newmat = B;
269: }
270: return(0);
271: }
273: PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)
274: {
275: Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data;
276: const PetscScalar *xv;
277: PetscScalar *yv;
278: PetscBLASInt N,m,ldax = 0,lday = 0,one = 1;
279: PetscErrorCode ierr;
282: MatDenseGetArrayRead(X,&xv);
283: MatDenseGetArray(Y,&yv);
284: PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);
285: PetscBLASIntCast(X->rmap->n,&m);
286: PetscBLASIntCast(x->lda,&ldax);
287: PetscBLASIntCast(y->lda,&lday);
288: if (ldax>m || lday>m) {
289: PetscInt j;
291: for (j=0; j<X->cmap->n; j++) {
292: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,xv+j*ldax,&one,yv+j*lday,&one));
293: }
294: } else {
295: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&alpha,xv,&one,yv,&one));
296: }
297: MatDenseRestoreArrayRead(X,&xv);
298: MatDenseRestoreArray(Y,&yv);
299: PetscLogFlops(PetscMax(2.0*N-1,0));
300: return(0);
301: }
303: static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
304: {
305: PetscLogDouble N = A->rmap->n*A->cmap->n;
308: info->block_size = 1.0;
309: info->nz_allocated = N;
310: info->nz_used = N;
311: info->nz_unneeded = 0;
312: info->assemblies = A->num_ass;
313: info->mallocs = 0;
314: info->memory = ((PetscObject)A)->mem;
315: info->fill_ratio_given = 0;
316: info->fill_ratio_needed = 0;
317: info->factor_mallocs = 0;
318: return(0);
319: }
321: PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
322: {
323: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
324: PetscScalar *v;
326: PetscBLASInt one = 1,j,nz,lda = 0;
329: MatDenseGetArray(A,&v);
330: PetscBLASIntCast(a->lda,&lda);
331: if (lda>A->rmap->n) {
332: PetscBLASIntCast(A->rmap->n,&nz);
333: for (j=0; j<A->cmap->n; j++) {
334: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v+j*lda,&one));
335: }
336: } else {
337: PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);
338: PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v,&one));
339: }
340: PetscLogFlops(nz);
341: MatDenseRestoreArray(A,&v);
342: return(0);
343: }
345: static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl)
346: {
347: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
348: PetscInt i,j,m = A->rmap->n,N = a->lda;
349: const PetscScalar *v;
350: PetscErrorCode ierr;
353: *fl = PETSC_FALSE;
354: if (A->rmap->n != A->cmap->n) return(0);
355: MatDenseGetArrayRead(A,&v);
356: for (i=0; i<m; i++) {
357: for (j=i; j<m; j++) {
358: if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) {
359: goto restore;
360: }
361: }
362: }
363: *fl = PETSC_TRUE;
364: restore:
365: MatDenseRestoreArrayRead(A,&v);
366: return(0);
367: }
369: static PetscErrorCode MatIsSymmetric_SeqDense(Mat A,PetscReal rtol,PetscBool *fl)
370: {
371: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
372: PetscInt i,j,m = A->rmap->n,N = a->lda;
373: const PetscScalar *v;
374: PetscErrorCode ierr;
377: *fl = PETSC_FALSE;
378: if (A->rmap->n != A->cmap->n) return(0);
379: MatDenseGetArrayRead(A,&v);
380: for (i=0; i<m; i++) {
381: for (j=i; j<m; j++) {
382: if (PetscAbsScalar(v[i+j*N] - v[j+i*N]) > rtol) {
383: goto restore;
384: }
385: }
386: }
387: *fl = PETSC_TRUE;
388: restore:
389: MatDenseRestoreArrayRead(A,&v);
390: return(0);
391: }
393: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
394: {
395: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
397: PetscInt lda = (PetscInt)mat->lda,j,m,nlda = lda;
400: PetscLayoutReference(A->rmap,&newi->rmap);
401: PetscLayoutReference(A->cmap,&newi->cmap);
402: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */
403: MatDenseSetLDA(newi,lda);
404: }
405: MatSeqDenseSetPreallocation(newi,NULL);
406: if (cpvalues == MAT_COPY_VALUES) {
407: const PetscScalar *av;
408: PetscScalar *v;
410: MatDenseGetArrayRead(A,&av);
411: MatDenseGetArray(newi,&v);
412: MatDenseGetLDA(newi,&nlda);
413: m = A->rmap->n;
414: if (lda>m || nlda>m) {
415: for (j=0; j<A->cmap->n; j++) {
416: PetscArraycpy(v+j*nlda,av+j*lda,m);
417: }
418: } else {
419: PetscArraycpy(v,av,A->rmap->n*A->cmap->n);
420: }
421: MatDenseRestoreArray(newi,&v);
422: MatDenseRestoreArrayRead(A,&av);
423: }
424: return(0);
425: }
427: PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
428: {
432: MatCreate(PetscObjectComm((PetscObject)A),newmat);
433: MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
434: MatSetType(*newmat,((PetscObject)A)->type_name);
435: MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);
436: return(0);
437: }
439: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
440: {
441: MatFactorInfo info;
445: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
446: (*fact->ops->lufactor)(fact,NULL,NULL,&info);
447: return(0);
448: }
450: static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy)
451: {
452: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
453: PetscErrorCode ierr;
454: const PetscScalar *x;
455: PetscScalar *y;
456: PetscBLASInt one = 1,info,m;
459: PetscBLASIntCast(A->rmap->n,&m);
460: VecGetArrayRead(xx,&x);
461: VecGetArray(yy,&y);
462: PetscArraycpy(y,x,A->rmap->n);
463: if (A->factortype == MAT_FACTOR_LU) {
464: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
465: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
466: PetscFPTrapPop();
467: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
468: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
469: if (A->spd) {
470: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
471: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
472: PetscFPTrapPop();
473: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
474: #if defined(PETSC_USE_COMPLEX)
475: } else if (A->hermitian) {
476: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
477: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
478: PetscFPTrapPop();
479: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
480: #endif
481: } else { /* symmetric case */
482: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
483: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
484: PetscFPTrapPop();
485: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
486: }
487: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
488: VecRestoreArrayRead(xx,&x);
489: VecRestoreArray(yy,&y);
490: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
491: return(0);
492: }
494: static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X)
495: {
496: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
497: PetscErrorCode ierr;
498: const PetscScalar *b;
499: PetscScalar *x;
500: PetscInt n;
501: PetscBLASInt nrhs,info,m;
504: PetscBLASIntCast(A->rmap->n,&m);
505: MatGetSize(B,NULL,&n);
506: PetscBLASIntCast(n,&nrhs);
507: MatDenseGetArrayRead(B,&b);
508: MatDenseGetArray(X,&x);
510: PetscArraycpy(x,b,m*nrhs);
512: if (A->factortype == MAT_FACTOR_LU) {
513: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
514: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
515: PetscFPTrapPop();
516: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve");
517: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
518: if (A->spd) {
519: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
520: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
521: PetscFPTrapPop();
522: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
523: #if defined(PETSC_USE_COMPLEX)
524: } else if (A->hermitian) {
525: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
526: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
527: PetscFPTrapPop();
528: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
529: #endif
530: } else { /* symmetric case */
531: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
532: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
533: PetscFPTrapPop();
534: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
535: }
536: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
538: MatDenseRestoreArrayRead(B,&b);
539: MatDenseRestoreArray(X,&x);
540: PetscLogFlops(nrhs*(2.0*m*m - m));
541: return(0);
542: }
544: static PetscErrorCode MatConjugate_SeqDense(Mat);
546: static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)
547: {
548: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
549: PetscErrorCode ierr;
550: const PetscScalar *x;
551: PetscScalar *y;
552: PetscBLASInt one = 1,info,m;
555: PetscBLASIntCast(A->rmap->n,&m);
556: VecGetArrayRead(xx,&x);
557: VecGetArray(yy,&y);
558: PetscArraycpy(y,x,A->rmap->n);
559: if (A->factortype == MAT_FACTOR_LU) {
560: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
561: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
562: PetscFPTrapPop();
563: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve");
564: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
565: if (A->spd) {
566: #if defined(PETSC_USE_COMPLEX)
567: MatConjugate_SeqDense(A);
568: #endif
569: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
570: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
571: PetscFPTrapPop();
572: #if defined(PETSC_USE_COMPLEX)
573: MatConjugate_SeqDense(A);
574: #endif
575: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve");
576: #if defined(PETSC_USE_COMPLEX)
577: } else if (A->hermitian) {
578: MatConjugate_SeqDense(A);
579: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
580: PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
581: PetscFPTrapPop();
582: MatConjugate_SeqDense(A);
583: #endif
584: } else { /* symmetric case */
585: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
586: PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
587: PetscFPTrapPop();
588: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve");
589: }
590: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
591: VecRestoreArrayRead(xx,&x);
592: VecRestoreArray(yy,&y);
593: PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);
594: return(0);
595: }
597: /* ---------------------------------------------------------------*/
598: /* COMMENT: I have chosen to hide row permutation in the pivots,
599: rather than put it in the Mat->row slot.*/
600: PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
601: {
602: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
604: PetscBLASInt n,m,info;
607: PetscBLASIntCast(A->cmap->n,&n);
608: PetscBLASIntCast(A->rmap->n,&m);
609: if (!mat->pivots) {
610: PetscMalloc1(A->rmap->n,&mat->pivots);
611: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
612: }
613: if (!A->rmap->n || !A->cmap->n) return(0);
614: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
615: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
616: PetscFPTrapPop();
618: if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization");
619: if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");
621: A->ops->solve = MatSolve_SeqDense;
622: A->ops->matsolve = MatMatSolve_SeqDense;
623: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
624: A->factortype = MAT_FACTOR_LU;
626: PetscFree(A->solvertype);
627: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
629: PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);
630: return(0);
631: }
633: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
634: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
635: {
636: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
638: PetscBLASInt info,n;
641: PetscBLASIntCast(A->cmap->n,&n);
642: if (!A->rmap->n || !A->cmap->n) return(0);
643: if (A->spd) {
644: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
645: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
646: PetscFPTrapPop();
647: #if defined(PETSC_USE_COMPLEX)
648: } else if (A->hermitian) {
649: if (!mat->pivots) {
650: PetscMalloc1(A->rmap->n,&mat->pivots);
651: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
652: }
653: if (!mat->fwork) {
654: PetscScalar dummy;
656: mat->lfwork = -1;
657: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
658: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
659: PetscFPTrapPop();
660: mat->lfwork = (PetscInt)PetscRealPart(dummy);
661: PetscMalloc1(mat->lfwork,&mat->fwork);
662: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
663: }
664: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
665: PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
666: PetscFPTrapPop();
667: #endif
668: } else { /* symmetric case */
669: if (!mat->pivots) {
670: PetscMalloc1(A->rmap->n,&mat->pivots);
671: PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));
672: }
673: if (!mat->fwork) {
674: PetscScalar dummy;
676: mat->lfwork = -1;
677: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
678: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
679: PetscFPTrapPop();
680: mat->lfwork = (PetscInt)PetscRealPart(dummy);
681: PetscMalloc1(mat->lfwork,&mat->fwork);
682: PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));
683: }
684: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
685: PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
686: PetscFPTrapPop();
687: }
688: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
690: A->ops->solve = MatSolve_SeqDense;
691: A->ops->matsolve = MatMatSolve_SeqDense;
692: A->ops->solvetranspose = MatSolveTranspose_SeqDense;
693: A->factortype = MAT_FACTOR_CHOLESKY;
695: PetscFree(A->solvertype);
696: PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);
698: PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);
699: return(0);
700: }
702: PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
703: {
705: MatFactorInfo info;
708: info.fill = 1.0;
710: MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);
711: (*fact->ops->choleskyfactor)(fact,NULL,&info);
712: return(0);
713: }
715: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
716: {
718: fact->assembled = PETSC_TRUE;
719: fact->preallocated = PETSC_TRUE;
720: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
721: fact->ops->solve = MatSolve_SeqDense;
722: fact->ops->matsolve = MatMatSolve_SeqDense;
723: fact->ops->solvetranspose = MatSolveTranspose_SeqDense;
724: return(0);
725: }
727: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
728: {
730: fact->preallocated = PETSC_TRUE;
731: fact->assembled = PETSC_TRUE;
732: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
733: fact->ops->solve = MatSolve_SeqDense;
734: fact->ops->matsolve = MatMatSolve_SeqDense;
735: fact->ops->solvetranspose = MatSolveTranspose_SeqDense;
736: return(0);
737: }
739: /* uses LAPACK */
740: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
741: {
745: MatCreate(PetscObjectComm((PetscObject)A),fact);
746: MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
747: MatSetType(*fact,MATDENSE);
748: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
749: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
750: (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
751: } else {
752: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
753: }
754: (*fact)->factortype = ftype;
756: PetscFree((*fact)->solvertype);
757: PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);
758: return(0);
759: }
761: /* ------------------------------------------------------------------*/
762: static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)
763: {
764: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
765: PetscScalar *x,*v = mat->v,zero = 0.0,xt;
766: const PetscScalar *b;
767: PetscErrorCode ierr;
768: PetscInt m = A->rmap->n,i;
769: PetscBLASInt o = 1,bm = 0;
772: #if defined(PETSC_HAVE_CUDA)
773: if (A->offloadmask == PETSC_OFFLOAD_GPU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented");
774: #endif
775: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
776: PetscBLASIntCast(m,&bm);
777: if (flag & SOR_ZERO_INITIAL_GUESS) {
778: /* this is a hack fix, should have another version without the second BLASdotu */
779: VecSet(xx,zero);
780: }
781: VecGetArray(xx,&x);
782: VecGetArrayRead(bb,&b);
783: its = its*lits;
784: if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
785: while (its--) {
786: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
787: for (i=0; i<m; i++) {
788: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
789: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
790: }
791: }
792: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
793: for (i=m-1; i>=0; i--) {
794: PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o));
795: x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift);
796: }
797: }
798: }
799: VecRestoreArrayRead(bb,&b);
800: VecRestoreArray(xx,&x);
801: return(0);
802: }
804: /* -----------------------------------------------------------------*/
805: PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)
806: {
807: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
808: const PetscScalar *v = mat->v,*x;
809: PetscScalar *y;
810: PetscErrorCode ierr;
811: PetscBLASInt m, n,_One=1;
812: PetscScalar _DOne=1.0,_DZero=0.0;
815: PetscBLASIntCast(A->rmap->n,&m);
816: PetscBLASIntCast(A->cmap->n,&n);
817: VecGetArrayRead(xx,&x);
818: VecGetArrayWrite(yy,&y);
819: if (!A->rmap->n || !A->cmap->n) {
820: PetscBLASInt i;
821: for (i=0; i<n; i++) y[i] = 0.0;
822: } else {
823: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One));
824: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);
825: }
826: VecRestoreArrayRead(xx,&x);
827: VecRestoreArrayWrite(yy,&y);
828: return(0);
829: }
831: PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy)
832: {
833: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
834: PetscScalar *y,_DOne=1.0,_DZero=0.0;
835: PetscErrorCode ierr;
836: PetscBLASInt m, n, _One=1;
837: const PetscScalar *v = mat->v,*x;
840: PetscBLASIntCast(A->rmap->n,&m);
841: PetscBLASIntCast(A->cmap->n,&n);
842: VecGetArrayRead(xx,&x);
843: VecGetArrayWrite(yy,&y);
844: if (!A->rmap->n || !A->cmap->n) {
845: PetscBLASInt i;
846: for (i=0; i<m; i++) y[i] = 0.0;
847: } else {
848: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One));
849: PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);
850: }
851: VecRestoreArrayRead(xx,&x);
852: VecRestoreArrayWrite(yy,&y);
853: return(0);
854: }
856: PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
857: {
858: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
859: const PetscScalar *v = mat->v,*x;
860: PetscScalar *y,_DOne=1.0;
861: PetscErrorCode ierr;
862: PetscBLASInt m, n, _One=1;
865: PetscBLASIntCast(A->rmap->n,&m);
866: PetscBLASIntCast(A->cmap->n,&n);
867: VecCopy(zz,yy);
868: if (!A->rmap->n || !A->cmap->n) return(0);
869: VecGetArrayRead(xx,&x);
870: VecGetArray(yy,&y);
871: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
872: VecRestoreArrayRead(xx,&x);
873: VecRestoreArray(yy,&y);
874: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
875: return(0);
876: }
878: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)
879: {
880: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
881: const PetscScalar *v = mat->v,*x;
882: PetscScalar *y;
883: PetscErrorCode ierr;
884: PetscBLASInt m, n, _One=1;
885: PetscScalar _DOne=1.0;
888: PetscBLASIntCast(A->rmap->n,&m);
889: PetscBLASIntCast(A->cmap->n,&n);
890: VecCopy(zz,yy);
891: if (!A->rmap->n || !A->cmap->n) return(0);
892: VecGetArrayRead(xx,&x);
893: VecGetArray(yy,&y);
894: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
895: VecRestoreArrayRead(xx,&x);
896: VecRestoreArray(yy,&y);
897: PetscLogFlops(2.0*A->rmap->n*A->cmap->n);
898: return(0);
899: }
901: /* -----------------------------------------------------------------*/
902: static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
903: {
904: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
906: PetscInt i;
909: *ncols = A->cmap->n;
910: if (cols) {
911: PetscMalloc1(A->cmap->n+1,cols);
912: for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
913: }
914: if (vals) {
915: const PetscScalar *v;
917: MatDenseGetArrayRead(A,&v);
918: PetscMalloc1(A->cmap->n+1,vals);
919: v += row;
920: for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
921: MatDenseRestoreArrayRead(A,&v);
922: }
923: return(0);
924: }
926: static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
927: {
931: if (cols) {PetscFree(*cols);}
932: if (vals) {PetscFree(*vals); }
933: return(0);
934: }
935: /* ----------------------------------------------------------------*/
936: static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)
937: {
938: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
939: PetscScalar *av;
940: PetscInt i,j,idx=0;
941: #if defined(PETSC_HAVE_CUDA)
942: PetscOffloadMask oldf;
943: #endif
944: PetscErrorCode ierr;
947: MatDenseGetArray(A,&av);
948: if (!mat->roworiented) {
949: if (addv == INSERT_VALUES) {
950: for (j=0; j<n; j++) {
951: if (indexn[j] < 0) {idx += m; continue;}
952: 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);
953: for (i=0; i<m; i++) {
954: if (indexm[i] < 0) {idx++; continue;}
955: 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);
956: av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
957: }
958: }
959: } else {
960: for (j=0; j<n; j++) {
961: if (indexn[j] < 0) {idx += m; continue;}
962: 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);
963: for (i=0; i<m; i++) {
964: if (indexm[i] < 0) {idx++; continue;}
965: 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);
966: av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
967: }
968: }
969: }
970: } else {
971: if (addv == INSERT_VALUES) {
972: for (i=0; i<m; i++) {
973: if (indexm[i] < 0) { idx += n; continue;}
974: 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);
975: for (j=0; j<n; j++) {
976: if (indexn[j] < 0) { idx++; continue;}
977: 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);
978: av[indexn[j]*mat->lda + indexm[i]] = v[idx++];
979: }
980: }
981: } else {
982: for (i=0; i<m; i++) {
983: if (indexm[i] < 0) { idx += n; continue;}
984: 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);
985: for (j=0; j<n; j++) {
986: if (indexn[j] < 0) { idx++; continue;}
987: 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);
988: av[indexn[j]*mat->lda + indexm[i]] += v[idx++];
989: }
990: }
991: }
992: }
993: /* hack to prevent unneeded copy to the GPU while returning the array */
994: #if defined(PETSC_HAVE_CUDA)
995: oldf = A->offloadmask;
996: A->offloadmask = PETSC_OFFLOAD_GPU;
997: #endif
998: MatDenseRestoreArray(A,&av);
999: #if defined(PETSC_HAVE_CUDA)
1000: A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1001: #endif
1002: return(0);
1003: }
1005: static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])
1006: {
1007: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1008: const PetscScalar *vv;
1009: PetscInt i,j;
1010: PetscErrorCode ierr;
1013: MatDenseGetArrayRead(A,&vv);
1014: /* row-oriented output */
1015: for (i=0; i<m; i++) {
1016: if (indexm[i] < 0) {v += n;continue;}
1017: 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);
1018: for (j=0; j<n; j++) {
1019: if (indexn[j] < 0) {v++; continue;}
1020: 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);
1021: *v++ = vv[indexn[j]*mat->lda + indexm[i]];
1022: }
1023: }
1024: MatDenseRestoreArrayRead(A,&vv);
1025: return(0);
1026: }
1028: /* -----------------------------------------------------------------*/
1030: PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer)
1031: {
1032: PetscErrorCode ierr;
1033: PetscBool skipHeader;
1034: PetscViewerFormat format;
1035: PetscInt header[4],M,N,m,lda,i,j,k;
1036: const PetscScalar *v;
1037: PetscScalar *vwork;
1040: PetscViewerSetUp(viewer);
1041: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
1042: PetscViewerGetFormat(viewer,&format);
1043: if (skipHeader) format = PETSC_VIEWER_NATIVE;
1045: MatGetSize(mat,&M,&N);
1047: /* write matrix header */
1048: header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N;
1049: header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N;
1050: if (!skipHeader) {PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);}
1052: MatGetLocalSize(mat,&m,NULL);
1053: if (format != PETSC_VIEWER_NATIVE) {
1054: PetscInt nnz = m*N, *iwork;
1055: /* store row lengths for each row */
1056: PetscMalloc1(nnz,&iwork);
1057: for (i=0; i<m; i++) iwork[i] = N;
1058: PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1059: /* store column indices (zero start index) */
1060: for (k=0, i=0; i<m; i++)
1061: for (j=0; j<N; j++, k++)
1062: iwork[k] = j;
1063: PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1064: PetscFree(iwork);
1065: }
1066: /* store matrix values as a dense matrix in row major order */
1067: PetscMalloc1(m*N,&vwork);
1068: MatDenseGetArrayRead(mat,&v);
1069: MatDenseGetLDA(mat,&lda);
1070: for (k=0, i=0; i<m; i++)
1071: for (j=0; j<N; j++, k++)
1072: vwork[k] = v[i+lda*j];
1073: MatDenseRestoreArrayRead(mat,&v);
1074: PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1075: PetscFree(vwork);
1076: return(0);
1077: }
1079: PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer)
1080: {
1082: PetscBool skipHeader;
1083: PetscInt header[4],M,N,m,nz,lda,i,j,k;
1084: PetscInt rows,cols;
1085: PetscScalar *v,*vwork;
1088: PetscViewerSetUp(viewer);
1089: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
1091: if (!skipHeader) {
1092: PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
1093: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file");
1094: M = header[1]; N = header[2];
1095: if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M);
1096: if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N);
1097: nz = header[3];
1098: if (nz != MATRIX_BINARY_FORMAT_DENSE && nz < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Unknown matrix format %D in file",nz);
1099: } else {
1100: MatGetSize(mat,&M,&N);
1101: 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");
1102: nz = MATRIX_BINARY_FORMAT_DENSE;
1103: }
1105: /* setup global sizes if not set */
1106: if (mat->rmap->N < 0) mat->rmap->N = M;
1107: if (mat->cmap->N < 0) mat->cmap->N = N;
1108: MatSetUp(mat);
1109: /* check if global sizes are correct */
1110: MatGetSize(mat,&rows,&cols);
1111: 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);
1113: MatGetSize(mat,NULL,&N);
1114: MatGetLocalSize(mat,&m,NULL);
1115: MatDenseGetArray(mat,&v);
1116: MatDenseGetLDA(mat,&lda);
1117: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */
1118: PetscInt nnz = m*N;
1119: /* read in matrix values */
1120: PetscMalloc1(nnz,&vwork);
1121: PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1122: /* store values in column major order */
1123: for (j=0; j<N; j++)
1124: for (i=0; i<m; i++)
1125: v[i+lda*j] = vwork[i*N+j];
1126: PetscFree(vwork);
1127: } else { /* matrix in file is sparse format */
1128: PetscInt nnz = 0, *rlens, *icols;
1129: /* read in row lengths */
1130: PetscMalloc1(m,&rlens);
1131: PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1132: for (i=0; i<m; i++) nnz += rlens[i];
1133: /* read in column indices and values */
1134: PetscMalloc2(nnz,&icols,nnz,&vwork);
1135: PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);
1136: PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);
1137: /* store values in column major order */
1138: for (k=0, i=0; i<m; i++)
1139: for (j=0; j<rlens[i]; j++, k++)
1140: v[i+lda*icols[k]] = vwork[k];
1141: PetscFree(rlens);
1142: PetscFree2(icols,vwork);
1143: }
1144: MatDenseRestoreArray(mat,&v);
1145: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
1146: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
1147: return(0);
1148: }
1150: PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1151: {
1152: PetscBool isbinary, ishdf5;
1158: /* force binary viewer to load .info file if it has not yet done so */
1159: PetscViewerSetUp(viewer);
1160: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1161: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);
1162: if (isbinary) {
1163: MatLoad_Dense_Binary(newMat,viewer);
1164: } else if (ishdf5) {
1165: #if defined(PETSC_HAVE_HDF5)
1166: MatLoad_Dense_HDF5(newMat,viewer);
1167: #else
1168: SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1169: #endif
1170: } else {
1171: 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);
1172: }
1173: return(0);
1174: }
1176: static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)
1177: {
1178: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1179: PetscErrorCode ierr;
1180: PetscInt i,j;
1181: const char *name;
1182: PetscScalar *v,*av;
1183: PetscViewerFormat format;
1184: #if defined(PETSC_USE_COMPLEX)
1185: PetscBool allreal = PETSC_TRUE;
1186: #endif
1189: MatDenseGetArrayRead(A,(const PetscScalar**)&av);
1190: PetscViewerGetFormat(viewer,&format);
1191: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1192: return(0); /* do nothing for now */
1193: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1194: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1195: for (i=0; i<A->rmap->n; i++) {
1196: v = av + i;
1197: PetscViewerASCIIPrintf(viewer,"row %D:",i);
1198: for (j=0; j<A->cmap->n; j++) {
1199: #if defined(PETSC_USE_COMPLEX)
1200: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1201: PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1202: } else if (PetscRealPart(*v)) {
1203: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));
1204: }
1205: #else
1206: if (*v) {
1207: PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);
1208: }
1209: #endif
1210: v += a->lda;
1211: }
1212: PetscViewerASCIIPrintf(viewer,"\n");
1213: }
1214: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1215: } else {
1216: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1217: #if defined(PETSC_USE_COMPLEX)
1218: /* determine if matrix has all real values */
1219: v = av;
1220: for (i=0; i<A->rmap->n*A->cmap->n; i++) {
1221: if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;}
1222: }
1223: #endif
1224: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1225: PetscObjectGetName((PetscObject)A,&name);
1226: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);
1227: PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);
1228: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
1229: }
1231: for (i=0; i<A->rmap->n; i++) {
1232: v = av + i;
1233: for (j=0; j<A->cmap->n; j++) {
1234: #if defined(PETSC_USE_COMPLEX)
1235: if (allreal) {
1236: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
1237: } else {
1238: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
1239: }
1240: #else
1241: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
1242: #endif
1243: v += a->lda;
1244: }
1245: PetscViewerASCIIPrintf(viewer,"\n");
1246: }
1247: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1248: PetscViewerASCIIPrintf(viewer,"];\n");
1249: }
1250: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1251: }
1252: MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);
1253: PetscViewerFlush(viewer);
1254: return(0);
1255: }
1257: #include <petscdraw.h>
1258: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa)
1259: {
1260: Mat A = (Mat) Aa;
1261: PetscErrorCode ierr;
1262: PetscInt m = A->rmap->n,n = A->cmap->n,i,j;
1263: int color = PETSC_DRAW_WHITE;
1264: const PetscScalar *v;
1265: PetscViewer viewer;
1266: PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1267: PetscViewerFormat format;
1270: PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1271: PetscViewerGetFormat(viewer,&format);
1272: PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
1274: /* Loop over matrix elements drawing boxes */
1275: MatDenseGetArrayRead(A,&v);
1276: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1277: PetscDrawCollectiveBegin(draw);
1278: /* Blue for negative and Red for positive */
1279: for (j = 0; j < n; j++) {
1280: x_l = j; x_r = x_l + 1.0;
1281: for (i = 0; i < m; i++) {
1282: y_l = m - i - 1.0;
1283: y_r = y_l + 1.0;
1284: if (PetscRealPart(v[j*m+i]) > 0.) color = PETSC_DRAW_RED;
1285: else if (PetscRealPart(v[j*m+i]) < 0.) color = PETSC_DRAW_BLUE;
1286: else continue;
1287: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1288: }
1289: }
1290: PetscDrawCollectiveEnd(draw);
1291: } else {
1292: /* use contour shading to indicate magnitude of values */
1293: /* first determine max of all nonzero values */
1294: PetscReal minv = 0.0, maxv = 0.0;
1295: PetscDraw popup;
1297: for (i=0; i < m*n; i++) {
1298: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1299: }
1300: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1301: PetscDrawGetPopup(draw,&popup);
1302: PetscDrawScalePopup(popup,minv,maxv);
1304: PetscDrawCollectiveBegin(draw);
1305: for (j=0; j<n; j++) {
1306: x_l = j;
1307: x_r = x_l + 1.0;
1308: for (i=0; i<m; i++) {
1309: y_l = m - i - 1.0;
1310: y_r = y_l + 1.0;
1311: color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv);
1312: PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);
1313: }
1314: }
1315: PetscDrawCollectiveEnd(draw);
1316: }
1317: MatDenseRestoreArrayRead(A,&v);
1318: return(0);
1319: }
1321: static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer)
1322: {
1323: PetscDraw draw;
1324: PetscBool isnull;
1325: PetscReal xr,yr,xl,yl,h,w;
1329: PetscViewerDrawGetDraw(viewer,0,&draw);
1330: PetscDrawIsNull(draw,&isnull);
1331: if (isnull) return(0);
1333: xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0;
1334: xr += w; yr += h; xl = -w; yl = -h;
1335: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1336: PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1337: PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);
1338: PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1339: PetscDrawSave(draw);
1340: return(0);
1341: }
1343: PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1344: {
1346: PetscBool iascii,isbinary,isdraw;
1349: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1350: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1351: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1353: if (iascii) {
1354: MatView_SeqDense_ASCII(A,viewer);
1355: } else if (isbinary) {
1356: MatView_Dense_Binary(A,viewer);
1357: } else if (isdraw) {
1358: MatView_SeqDense_Draw(A,viewer);
1359: }
1360: return(0);
1361: }
1363: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar *array)
1364: {
1365: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1368: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1369: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1370: if (a->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreArray() first");
1371: a->unplacedarray = a->v;
1372: a->unplaced_user_alloc = a->user_alloc;
1373: a->v = (PetscScalar*) array;
1374: a->user_alloc = PETSC_TRUE;
1375: #if defined(PETSC_HAVE_CUDA)
1376: A->offloadmask = PETSC_OFFLOAD_CPU;
1377: #endif
1378: return(0);
1379: }
1381: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1382: {
1383: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1386: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1387: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1388: a->v = a->unplacedarray;
1389: a->user_alloc = a->unplaced_user_alloc;
1390: a->unplacedarray = NULL;
1391: #if defined(PETSC_HAVE_CUDA)
1392: A->offloadmask = PETSC_OFFLOAD_CPU;
1393: #endif
1394: return(0);
1395: }
1397: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A,const PetscScalar *array)
1398: {
1399: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1403: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1404: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1405: if (!a->user_alloc) { PetscFree(a->v); }
1406: a->v = (PetscScalar*) array;
1407: a->user_alloc = PETSC_FALSE;
1408: #if defined(PETSC_HAVE_CUDA)
1409: A->offloadmask = PETSC_OFFLOAD_CPU;
1410: #endif
1411: return(0);
1412: }
1414: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1415: {
1416: Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
1420: #if defined(PETSC_USE_LOG)
1421: PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1422: #endif
1423: PetscFree(l->pivots);
1424: PetscFree(l->fwork);
1425: MatDestroy(&l->ptapwork);
1426: if (!l->user_alloc) {PetscFree(l->v);}
1427: if (!l->unplaced_user_alloc) {PetscFree(l->unplacedarray);}
1428: if (l->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
1429: if (l->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1430: VecDestroy(&l->cvec);
1431: MatDestroy(&l->cmat);
1432: PetscFree(mat->data);
1434: PetscObjectChangeTypeName((PetscObject)mat,NULL);
1435: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);
1436: PetscObjectComposeFunction((PetscObject)mat,"MatDenseSetLDA_C",NULL);
1437: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);
1438: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);
1439: PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);
1440: PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);
1441: PetscObjectComposeFunction((PetscObject)mat,"MatDenseReplaceArray_C",NULL);
1442: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);
1443: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);
1444: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayWrite_C",NULL);
1445: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayWrite_C",NULL);
1446: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);
1447: #if defined(PETSC_HAVE_ELEMENTAL)
1448: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);
1449: #endif
1450: #if defined(PETSC_HAVE_SCALAPACK)
1451: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_scalapack_C",NULL);
1452: #endif
1453: #if defined(PETSC_HAVE_CUDA)
1454: PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);
1455: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL);
1456: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL);
1457: #endif
1458: PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);
1459: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL);
1460: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL);
1461: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL);
1462: PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL);
1464: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);
1465: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);
1466: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVec_C",NULL);
1467: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVec_C",NULL);
1468: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecRead_C",NULL);
1469: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecRead_C",NULL);
1470: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecWrite_C",NULL);
1471: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecWrite_C",NULL);
1472: PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetSubMatrix_C",NULL);
1473: PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreSubMatrix_C",NULL);
1474: return(0);
1475: }
1477: static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1478: {
1479: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1481: PetscInt k,j,m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1482: PetscScalar *v,tmp;
1485: if (reuse == MAT_INPLACE_MATRIX) {
1486: if (m == n) { /* in place transpose */
1487: MatDenseGetArray(A,&v);
1488: for (j=0; j<m; j++) {
1489: for (k=0; k<j; k++) {
1490: tmp = v[j + k*M];
1491: v[j + k*M] = v[k + j*M];
1492: v[k + j*M] = tmp;
1493: }
1494: }
1495: MatDenseRestoreArray(A,&v);
1496: } else { /* reuse memory, temporary allocates new memory */
1497: PetscScalar *v2;
1498: PetscLayout tmplayout;
1500: PetscMalloc1((size_t)m*n,&v2);
1501: MatDenseGetArray(A,&v);
1502: for (j=0; j<n; j++) {
1503: for (k=0; k<m; k++) v2[j + (size_t)k*n] = v[k + (size_t)j*M];
1504: }
1505: PetscArraycpy(v,v2,(size_t)m*n);
1506: PetscFree(v2);
1507: MatDenseRestoreArray(A,&v);
1508: /* cleanup size dependent quantities */
1509: VecDestroy(&mat->cvec);
1510: MatDestroy(&mat->cmat);
1511: PetscFree(mat->pivots);
1512: PetscFree(mat->fwork);
1513: MatDestroy(&mat->ptapwork);
1514: /* swap row/col layouts */
1515: mat->lda = n;
1516: tmplayout = A->rmap;
1517: A->rmap = A->cmap;
1518: A->cmap = tmplayout;
1519: }
1520: } else { /* out-of-place transpose */
1521: Mat tmat;
1522: Mat_SeqDense *tmatd;
1523: PetscScalar *v2;
1524: PetscInt M2;
1526: if (reuse == MAT_INITIAL_MATRIX) {
1527: MatCreate(PetscObjectComm((PetscObject)A),&tmat);
1528: MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);
1529: MatSetType(tmat,((PetscObject)A)->type_name);
1530: MatSeqDenseSetPreallocation(tmat,NULL);
1531: } else tmat = *matout;
1533: MatDenseGetArrayRead(A,(const PetscScalar**)&v);
1534: MatDenseGetArray(tmat,&v2);
1535: tmatd = (Mat_SeqDense*)tmat->data;
1536: M2 = tmatd->lda;
1537: for (j=0; j<n; j++) {
1538: for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M];
1539: }
1540: MatDenseRestoreArray(tmat,&v2);
1541: MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);
1542: MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);
1543: MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);
1544: *matout = tmat;
1545: }
1546: return(0);
1547: }
1549: static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg)
1550: {
1551: Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data;
1552: Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data;
1553: PetscInt i;
1554: const PetscScalar *v1,*v2;
1555: PetscErrorCode ierr;
1558: if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; return(0);}
1559: if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; return(0);}
1560: MatDenseGetArrayRead(A1,&v1);
1561: MatDenseGetArrayRead(A2,&v2);
1562: for (i=0; i<A1->cmap->n; i++) {
1563: PetscArraycmp(v1,v2,A1->rmap->n,flg);
1564: if (*flg == PETSC_FALSE) return(0);
1565: v1 += mat1->lda;
1566: v2 += mat2->lda;
1567: }
1568: MatDenseRestoreArrayRead(A1,&v1);
1569: MatDenseRestoreArrayRead(A2,&v2);
1570: *flg = PETSC_TRUE;
1571: return(0);
1572: }
1574: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v)
1575: {
1576: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1577: PetscInt i,n,len;
1578: PetscScalar *x;
1579: const PetscScalar *vv;
1580: PetscErrorCode ierr;
1583: VecGetSize(v,&n);
1584: VecGetArray(v,&x);
1585: len = PetscMin(A->rmap->n,A->cmap->n);
1586: MatDenseGetArrayRead(A,&vv);
1587: if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec");
1588: for (i=0; i<len; i++) {
1589: x[i] = vv[i*mat->lda + i];
1590: }
1591: MatDenseRestoreArrayRead(A,&vv);
1592: VecRestoreArray(v,&x);
1593: return(0);
1594: }
1596: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)
1597: {
1598: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1599: const PetscScalar *l,*r;
1600: PetscScalar x,*v,*vv;
1601: PetscErrorCode ierr;
1602: PetscInt i,j,m = A->rmap->n,n = A->cmap->n;
1605: MatDenseGetArray(A,&vv);
1606: if (ll) {
1607: VecGetSize(ll,&m);
1608: VecGetArrayRead(ll,&l);
1609: if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size");
1610: for (i=0; i<m; i++) {
1611: x = l[i];
1612: v = vv + i;
1613: for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;}
1614: }
1615: VecRestoreArrayRead(ll,&l);
1616: PetscLogFlops(1.0*n*m);
1617: }
1618: if (rr) {
1619: VecGetSize(rr,&n);
1620: VecGetArrayRead(rr,&r);
1621: if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size");
1622: for (i=0; i<n; i++) {
1623: x = r[i];
1624: v = vv + i*mat->lda;
1625: for (j=0; j<m; j++) (*v++) *= x;
1626: }
1627: VecRestoreArrayRead(rr,&r);
1628: PetscLogFlops(1.0*n*m);
1629: }
1630: MatDenseRestoreArray(A,&vv);
1631: return(0);
1632: }
1634: PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm)
1635: {
1636: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1637: PetscScalar *v,*vv;
1638: PetscReal sum = 0.0;
1639: PetscInt lda =mat->lda,m=A->rmap->n,i,j;
1640: PetscErrorCode ierr;
1643: MatDenseGetArrayRead(A,(const PetscScalar**)&vv);
1644: v = vv;
1645: if (type == NORM_FROBENIUS) {
1646: if (lda>m) {
1647: for (j=0; j<A->cmap->n; j++) {
1648: v = vv+j*lda;
1649: for (i=0; i<m; i++) {
1650: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1651: }
1652: }
1653: } else {
1654: #if defined(PETSC_USE_REAL___FP16)
1655: PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n;
1656: *nrm = BLASnrm2_(&cnt,v,&one);
1657: }
1658: #else
1659: for (i=0; i<A->cmap->n*A->rmap->n; i++) {
1660: sum += PetscRealPart(PetscConj(*v)*(*v)); v++;
1661: }
1662: }
1663: *nrm = PetscSqrtReal(sum);
1664: #endif
1665: PetscLogFlops(2.0*A->cmap->n*A->rmap->n);
1666: } else if (type == NORM_1) {
1667: *nrm = 0.0;
1668: for (j=0; j<A->cmap->n; j++) {
1669: v = vv + j*mat->lda;
1670: sum = 0.0;
1671: for (i=0; i<A->rmap->n; i++) {
1672: sum += PetscAbsScalar(*v); v++;
1673: }
1674: if (sum > *nrm) *nrm = sum;
1675: }
1676: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1677: } else if (type == NORM_INFINITY) {
1678: *nrm = 0.0;
1679: for (j=0; j<A->rmap->n; j++) {
1680: v = vv + j;
1681: sum = 0.0;
1682: for (i=0; i<A->cmap->n; i++) {
1683: sum += PetscAbsScalar(*v); v += mat->lda;
1684: }
1685: if (sum > *nrm) *nrm = sum;
1686: }
1687: PetscLogFlops(1.0*A->cmap->n*A->rmap->n);
1688: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1689: MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);
1690: return(0);
1691: }
1693: static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1694: {
1695: Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
1699: switch (op) {
1700: case MAT_ROW_ORIENTED:
1701: aij->roworiented = flg;
1702: break;
1703: case MAT_NEW_NONZERO_LOCATIONS:
1704: case MAT_NEW_NONZERO_LOCATION_ERR:
1705: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1706: case MAT_NEW_DIAGONALS:
1707: case MAT_KEEP_NONZERO_PATTERN:
1708: case MAT_IGNORE_OFF_PROC_ENTRIES:
1709: case MAT_USE_HASH_TABLE:
1710: case MAT_IGNORE_ZERO_ENTRIES:
1711: case MAT_IGNORE_LOWER_TRIANGULAR:
1712: case MAT_SORTED_FULL:
1713: PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1714: break;
1715: case MAT_SPD:
1716: case MAT_SYMMETRIC:
1717: case MAT_STRUCTURALLY_SYMMETRIC:
1718: case MAT_HERMITIAN:
1719: case MAT_SYMMETRY_ETERNAL:
1720: /* These options are handled directly by MatSetOption() */
1721: break;
1722: default:
1723: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]);
1724: }
1725: return(0);
1726: }
1728: static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1729: {
1730: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1732: PetscInt lda=l->lda,m=A->rmap->n,j;
1733: PetscScalar *v;
1736: MatDenseGetArray(A,&v);
1737: if (lda>m) {
1738: for (j=0; j<A->cmap->n; j++) {
1739: PetscArrayzero(v+j*lda,m);
1740: }
1741: } else {
1742: PetscArrayzero(v,A->rmap->n*A->cmap->n);
1743: }
1744: MatDenseRestoreArray(A,&v);
1745: return(0);
1746: }
1748: static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1749: {
1750: PetscErrorCode ierr;
1751: Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1752: PetscInt m = l->lda, n = A->cmap->n, i,j;
1753: PetscScalar *slot,*bb,*v;
1754: const PetscScalar *xx;
1757: if (PetscDefined(USE_DEBUG)) {
1758: for (i=0; i<N; i++) {
1759: if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed");
1760: 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);
1761: }
1762: }
1763: if (!N) return(0);
1765: /* fix right hand side if needed */
1766: if (x && b) {
1767: VecGetArrayRead(x,&xx);
1768: VecGetArray(b,&bb);
1769: for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1770: VecRestoreArrayRead(x,&xx);
1771: VecRestoreArray(b,&bb);
1772: }
1774: MatDenseGetArray(A,&v);
1775: for (i=0; i<N; i++) {
1776: slot = v + rows[i];
1777: for (j=0; j<n; j++) { *slot = 0.0; slot += m;}
1778: }
1779: if (diag != 0.0) {
1780: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
1781: for (i=0; i<N; i++) {
1782: slot = v + (m+1)*rows[i];
1783: *slot = diag;
1784: }
1785: }
1786: MatDenseRestoreArray(A,&v);
1787: return(0);
1788: }
1790: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda)
1791: {
1792: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1795: *lda = mat->lda;
1796: return(0);
1797: }
1799: PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array)
1800: {
1801: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1804: if (mat->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1805: *array = mat->v;
1806: return(0);
1807: }
1809: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array)
1810: {
1812: *array = NULL;
1813: return(0);
1814: }
1816: /*@C
1817: MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray()
1819: Not collective
1821: Input Parameter:
1822: . mat - a MATSEQDENSE or MATMPIDENSE matrix
1824: Output Parameter:
1825: . lda - the leading dimension
1827: Level: intermediate
1829: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseSetLDA()
1830: @*/
1831: PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda)
1832: {
1838: PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));
1839: return(0);
1840: }
1842: /*@C
1843: MatDenseSetLDA - Sets the leading dimension of the array used by the dense matrix
1845: Not collective
1847: Input Parameter:
1848: + mat - a MATSEQDENSE or MATMPIDENSE matrix
1849: - lda - the leading dimension
1851: Level: intermediate
1853: .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetLDA()
1854: @*/
1855: PetscErrorCode MatDenseSetLDA(Mat A,PetscInt lda)
1856: {
1861: PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda));
1862: return(0);
1863: }
1865: /*@C
1866: MatDenseGetArray - gives read-write access to the array where the data for a dense matrix is stored
1868: Logically Collective on Mat
1870: Input Parameter:
1871: . mat - a dense matrix
1873: Output Parameter:
1874: . array - pointer to the data
1876: Level: intermediate
1878: .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
1879: @*/
1880: PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array)
1881: {
1887: PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));
1888: return(0);
1889: }
1891: /*@C
1892: MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
1894: Logically Collective on Mat
1896: Input Parameters:
1897: + mat - a dense matrix
1898: - array - pointer to the data
1900: Level: intermediate
1902: .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
1903: @*/
1904: PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array)
1905: {
1911: PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));
1912: PetscObjectStateIncrease((PetscObject)A);
1913: #if defined(PETSC_HAVE_CUDA)
1914: A->offloadmask = PETSC_OFFLOAD_CPU;
1915: #endif
1916: return(0);
1917: }
1919: /*@C
1920: MatDenseGetArrayRead - gives read-only access to the array where the data for a dense matrix is stored
1922: Not Collective
1924: Input Parameter:
1925: . mat - a dense matrix
1927: Output Parameter:
1928: . array - pointer to the data
1930: Level: intermediate
1932: .seealso: MatDenseRestoreArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
1933: @*/
1934: PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array)
1935: {
1941: PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));
1942: return(0);
1943: }
1945: /*@C
1946: MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayRead()
1948: Not Collective
1950: Input Parameters:
1951: + mat - a dense matrix
1952: - array - pointer to the data
1954: Level: intermediate
1956: .seealso: MatDenseGetArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
1957: @*/
1958: PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
1959: {
1965: PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));
1966: return(0);
1967: }
1969: /*@C
1970: MatDenseGetArrayWrite - gives write-only access to the array where the data for a dense matrix is stored
1972: Not Collective
1974: Input Parameter:
1975: . mat - a dense matrix
1977: Output Parameter:
1978: . array - pointer to the data
1980: Level: intermediate
1982: .seealso: MatDenseRestoreArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
1983: @*/
1984: PetscErrorCode MatDenseGetArrayWrite(Mat A,PetscScalar **array)
1985: {
1991: PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array));
1992: return(0);
1993: }
1995: /*@C
1996: MatDenseRestoreArrayWrite - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayWrite()
1998: Not Collective
2000: Input Parameters:
2001: + mat - a dense matrix
2002: - array - pointer to the data
2004: Level: intermediate
2006: .seealso: MatDenseGetArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
2007: @*/
2008: PetscErrorCode MatDenseRestoreArrayWrite(Mat A,PetscScalar **array)
2009: {
2015: PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array));
2016: PetscObjectStateIncrease((PetscObject)A);
2017: #if defined(PETSC_HAVE_CUDA)
2018: A->offloadmask = PETSC_OFFLOAD_CPU;
2019: #endif
2020: return(0);
2021: }
2023: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B)
2024: {
2025: Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
2027: PetscInt i,j,nrows,ncols,blda;
2028: const PetscInt *irow,*icol;
2029: PetscScalar *av,*bv,*v = mat->v;
2030: Mat newmat;
2033: ISGetIndices(isrow,&irow);
2034: ISGetIndices(iscol,&icol);
2035: ISGetLocalSize(isrow,&nrows);
2036: ISGetLocalSize(iscol,&ncols);
2038: /* Check submatrixcall */
2039: if (scall == MAT_REUSE_MATRIX) {
2040: PetscInt n_cols,n_rows;
2041: MatGetSize(*B,&n_rows,&n_cols);
2042: if (n_rows != nrows || n_cols != ncols) {
2043: /* resize the result matrix to match number of requested rows/columns */
2044: MatSetSizes(*B,nrows,ncols,nrows,ncols);
2045: }
2046: newmat = *B;
2047: } else {
2048: /* Create and fill new matrix */
2049: MatCreate(PetscObjectComm((PetscObject)A),&newmat);
2050: MatSetSizes(newmat,nrows,ncols,nrows,ncols);
2051: MatSetType(newmat,((PetscObject)A)->type_name);
2052: MatSeqDenseSetPreallocation(newmat,NULL);
2053: }
2055: /* Now extract the data pointers and do the copy,column at a time */
2056: MatDenseGetArray(newmat,&bv);
2057: MatDenseGetLDA(newmat,&blda);
2058: for (i=0; i<ncols; i++) {
2059: av = v + mat->lda*icol[i];
2060: for (j=0; j<nrows; j++) bv[j] = av[irow[j]];
2061: bv += blda;
2062: }
2063: MatDenseRestoreArray(newmat,&bv);
2065: /* Assemble the matrices so that the correct flags are set */
2066: MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
2067: MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
2069: /* Free work space */
2070: ISRestoreIndices(isrow,&irow);
2071: ISRestoreIndices(iscol,&icol);
2072: *B = newmat;
2073: return(0);
2074: }
2076: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
2077: {
2079: PetscInt i;
2082: if (scall == MAT_INITIAL_MATRIX) {
2083: PetscCalloc1(n+1,B);
2084: }
2086: for (i=0; i<n; i++) {
2087: MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);
2088: }
2089: return(0);
2090: }
2092: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
2093: {
2095: return(0);
2096: }
2098: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
2099: {
2101: return(0);
2102: }
2104: static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str)
2105: {
2106: Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data;
2107: PetscErrorCode ierr;
2108: const PetscScalar *va;
2109: PetscScalar *vb;
2110: PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j;
2113: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2114: if (A->ops->copy != B->ops->copy) {
2115: MatCopy_Basic(A,B,str);
2116: return(0);
2117: }
2118: if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)");
2119: MatDenseGetArrayRead(A,&va);
2120: MatDenseGetArray(B,&vb);
2121: if (lda1>m || lda2>m) {
2122: for (j=0; j<n; j++) {
2123: PetscArraycpy(vb+j*lda2,va+j*lda1,m);
2124: }
2125: } else {
2126: PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);
2127: }
2128: MatDenseRestoreArray(B,&vb);
2129: MatDenseRestoreArrayRead(A,&va);
2130: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2131: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2132: return(0);
2133: }
2135: static PetscErrorCode MatSetUp_SeqDense(Mat A)
2136: {
2140: PetscLayoutSetUp(A->rmap);
2141: PetscLayoutSetUp(A->cmap);
2142: if (!A->preallocated) {
2143: MatSeqDenseSetPreallocation(A,NULL);
2144: }
2145: return(0);
2146: }
2148: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2149: {
2150: PetscInt i,nz = A->rmap->n*A->cmap->n;
2151: PetscScalar *aa;
2155: MatDenseGetArray(A,&aa);
2156: for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
2157: MatDenseRestoreArray(A,&aa);
2158: return(0);
2159: }
2161: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2162: {
2163: PetscInt i,nz = A->rmap->n*A->cmap->n;
2164: PetscScalar *aa;
2168: MatDenseGetArray(A,&aa);
2169: for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
2170: MatDenseRestoreArray(A,&aa);
2171: return(0);
2172: }
2174: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2175: {
2176: PetscInt i,nz = A->rmap->n*A->cmap->n;
2177: PetscScalar *aa;
2181: MatDenseGetArray(A,&aa);
2182: for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2183: MatDenseRestoreArray(A,&aa);
2184: return(0);
2185: }
2187: /* ----------------------------------------------------------------*/
2188: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2189: {
2191: PetscInt m=A->rmap->n,n=B->cmap->n;
2192: PetscBool cisdense;
2195: MatSetSizes(C,m,n,m,n);
2196: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2197: if (!cisdense) {
2198: PetscBool flg;
2200: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2201: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2202: }
2203: MatSetUp(C);
2204: return(0);
2205: }
2207: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2208: {
2209: Mat_SeqDense *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data;
2210: PetscBLASInt m,n,k;
2211: const PetscScalar *av,*bv;
2212: PetscScalar *cv;
2213: PetscScalar _DOne=1.0,_DZero=0.0;
2214: PetscErrorCode ierr;
2217: PetscBLASIntCast(C->rmap->n,&m);
2218: PetscBLASIntCast(C->cmap->n,&n);
2219: PetscBLASIntCast(A->cmap->n,&k);
2220: if (!m || !n || !k) return(0);
2221: MatDenseGetArrayRead(A,&av);
2222: MatDenseGetArrayRead(B,&bv);
2223: MatDenseGetArrayWrite(C,&cv);
2224: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2225: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2226: MatDenseRestoreArrayRead(A,&av);
2227: MatDenseRestoreArrayRead(B,&bv);
2228: MatDenseRestoreArrayWrite(C,&cv);
2229: return(0);
2230: }
2232: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2233: {
2235: PetscInt m=A->rmap->n,n=B->rmap->n;
2236: PetscBool cisdense;
2239: MatSetSizes(C,m,n,m,n);
2240: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2241: if (!cisdense) {
2242: PetscBool flg;
2244: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2245: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2246: }
2247: MatSetUp(C);
2248: return(0);
2249: }
2251: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2252: {
2253: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2254: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2255: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2256: const PetscScalar *av,*bv;
2257: PetscScalar *cv;
2258: PetscBLASInt m,n,k;
2259: PetscScalar _DOne=1.0,_DZero=0.0;
2260: PetscErrorCode ierr;
2263: PetscBLASIntCast(C->rmap->n,&m);
2264: PetscBLASIntCast(C->cmap->n,&n);
2265: PetscBLASIntCast(A->cmap->n,&k);
2266: if (!m || !n || !k) return(0);
2267: MatDenseGetArrayRead(A,&av);
2268: MatDenseGetArrayRead(B,&bv);
2269: MatDenseGetArrayWrite(C,&cv);
2270: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2271: MatDenseRestoreArrayRead(A,&av);
2272: MatDenseRestoreArrayRead(B,&bv);
2273: MatDenseRestoreArrayWrite(C,&cv);
2274: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2275: return(0);
2276: }
2278: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2279: {
2281: PetscInt m=A->cmap->n,n=B->cmap->n;
2282: PetscBool cisdense;
2285: MatSetSizes(C,m,n,m,n);
2286: PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");
2287: if (!cisdense) {
2288: PetscBool flg;
2290: PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);
2291: MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);
2292: }
2293: MatSetUp(C);
2294: return(0);
2295: }
2297: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)
2298: {
2299: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2300: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2301: Mat_SeqDense *c = (Mat_SeqDense*)C->data;
2302: const PetscScalar *av,*bv;
2303: PetscScalar *cv;
2304: PetscBLASInt m,n,k;
2305: PetscScalar _DOne=1.0,_DZero=0.0;
2306: PetscErrorCode ierr;
2309: PetscBLASIntCast(C->rmap->n,&m);
2310: PetscBLASIntCast(C->cmap->n,&n);
2311: PetscBLASIntCast(A->rmap->n,&k);
2312: if (!m || !n || !k) return(0);
2313: MatDenseGetArrayRead(A,&av);
2314: MatDenseGetArrayRead(B,&bv);
2315: MatDenseGetArrayWrite(C,&cv);
2316: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2317: MatDenseRestoreArrayRead(A,&av);
2318: MatDenseRestoreArrayRead(B,&bv);
2319: MatDenseRestoreArrayWrite(C,&cv);
2320: PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));
2321: return(0);
2322: }
2324: /* ----------------------------------------------- */
2325: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2326: {
2328: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2329: C->ops->productsymbolic = MatProductSymbolic_AB;
2330: return(0);
2331: }
2333: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2334: {
2336: C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2337: C->ops->productsymbolic = MatProductSymbolic_AtB;
2338: return(0);
2339: }
2341: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2342: {
2344: C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2345: C->ops->productsymbolic = MatProductSymbolic_ABt;
2346: return(0);
2347: }
2349: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2350: {
2352: Mat_Product *product = C->product;
2355: switch (product->type) {
2356: case MATPRODUCT_AB:
2357: MatProductSetFromOptions_SeqDense_AB(C);
2358: break;
2359: case MATPRODUCT_AtB:
2360: MatProductSetFromOptions_SeqDense_AtB(C);
2361: break;
2362: case MATPRODUCT_ABt:
2363: MatProductSetFromOptions_SeqDense_ABt(C);
2364: break;
2365: default:
2366: break;
2367: }
2368: return(0);
2369: }
2370: /* ----------------------------------------------- */
2372: static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])
2373: {
2374: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2375: PetscErrorCode ierr;
2376: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2377: PetscScalar *x;
2378: const PetscScalar *aa;
2381: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2382: VecGetArray(v,&x);
2383: VecGetLocalSize(v,&p);
2384: MatDenseGetArrayRead(A,&aa);
2385: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2386: for (i=0; i<m; i++) {
2387: x[i] = aa[i]; if (idx) idx[i] = 0;
2388: for (j=1; j<n; j++) {
2389: if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2390: }
2391: }
2392: MatDenseRestoreArrayRead(A,&aa);
2393: VecRestoreArray(v,&x);
2394: return(0);
2395: }
2397: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])
2398: {
2399: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2400: PetscErrorCode ierr;
2401: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2402: PetscScalar *x;
2403: PetscReal atmp;
2404: const PetscScalar *aa;
2407: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2408: VecGetArray(v,&x);
2409: VecGetLocalSize(v,&p);
2410: MatDenseGetArrayRead(A,&aa);
2411: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2412: for (i=0; i<m; i++) {
2413: x[i] = PetscAbsScalar(aa[i]);
2414: for (j=1; j<n; j++) {
2415: atmp = PetscAbsScalar(aa[i+a->lda*j]);
2416: if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;}
2417: }
2418: }
2419: MatDenseRestoreArrayRead(A,&aa);
2420: VecRestoreArray(v,&x);
2421: return(0);
2422: }
2424: static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])
2425: {
2426: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2427: PetscErrorCode ierr;
2428: PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p;
2429: PetscScalar *x;
2430: const PetscScalar *aa;
2433: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2434: MatDenseGetArrayRead(A,&aa);
2435: VecGetArray(v,&x);
2436: VecGetLocalSize(v,&p);
2437: if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2438: for (i=0; i<m; i++) {
2439: x[i] = aa[i]; if (idx) idx[i] = 0;
2440: for (j=1; j<n; j++) {
2441: if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;}
2442: }
2443: }
2444: VecRestoreArray(v,&x);
2445: MatDenseRestoreArrayRead(A,&aa);
2446: return(0);
2447: }
2449: PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)
2450: {
2451: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2452: PetscErrorCode ierr;
2453: PetscScalar *x;
2454: const PetscScalar *aa;
2457: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2458: MatDenseGetArrayRead(A,&aa);
2459: VecGetArray(v,&x);
2460: PetscArraycpy(x,aa+col*a->lda,A->rmap->n);
2461: VecRestoreArray(v,&x);
2462: MatDenseRestoreArrayRead(A,&aa);
2463: return(0);
2464: }
2466: PETSC_INTERN PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms)
2467: {
2468: PetscErrorCode ierr;
2469: PetscInt i,j,m,n;
2470: const PetscScalar *a;
2473: MatGetSize(A,&m,&n);
2474: PetscArrayzero(norms,n);
2475: MatDenseGetArrayRead(A,&a);
2476: if (type == NORM_2) {
2477: for (i=0; i<n; i++) {
2478: for (j=0; j<m; j++) {
2479: norms[i] += PetscAbsScalar(a[j]*a[j]);
2480: }
2481: a += m;
2482: }
2483: } else if (type == NORM_1) {
2484: for (i=0; i<n; i++) {
2485: for (j=0; j<m; j++) {
2486: norms[i] += PetscAbsScalar(a[j]);
2487: }
2488: a += m;
2489: }
2490: } else if (type == NORM_INFINITY) {
2491: for (i=0; i<n; i++) {
2492: for (j=0; j<m; j++) {
2493: norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]);
2494: }
2495: a += m;
2496: }
2497: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType");
2498: MatDenseRestoreArrayRead(A,&a);
2499: if (type == NORM_2) {
2500: for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2501: }
2502: return(0);
2503: }
2505: static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2506: {
2508: PetscScalar *a;
2509: PetscInt lda,m,n,i,j;
2512: MatGetSize(x,&m,&n);
2513: MatDenseGetLDA(x,&lda);
2514: MatDenseGetArray(x,&a);
2515: for (j=0; j<n; j++) {
2516: for (i=0; i<m; i++) {
2517: PetscRandomGetValue(rctx,a+j*lda+i);
2518: }
2519: }
2520: MatDenseRestoreArray(x,&a);
2521: return(0);
2522: }
2524: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d)
2525: {
2527: *missing = PETSC_FALSE;
2528: return(0);
2529: }
2531: /* vals is not const */
2532: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2533: {
2535: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2536: PetscScalar *v;
2539: if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2540: MatDenseGetArray(A,&v);
2541: *vals = v+col*a->lda;
2542: MatDenseRestoreArray(A,&v);
2543: return(0);
2544: }
2546: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2547: {
2549: *vals = NULL; /* user cannot accidently use the array later */
2550: return(0);
2551: }
2553: /* -------------------------------------------------------------------*/
2554: static struct _MatOps MatOps_Values = { MatSetValues_SeqDense,
2555: MatGetRow_SeqDense,
2556: MatRestoreRow_SeqDense,
2557: MatMult_SeqDense,
2558: /* 4*/ MatMultAdd_SeqDense,
2559: MatMultTranspose_SeqDense,
2560: MatMultTransposeAdd_SeqDense,
2561: NULL,
2562: NULL,
2563: NULL,
2564: /* 10*/ NULL,
2565: MatLUFactor_SeqDense,
2566: MatCholeskyFactor_SeqDense,
2567: MatSOR_SeqDense,
2568: MatTranspose_SeqDense,
2569: /* 15*/ MatGetInfo_SeqDense,
2570: MatEqual_SeqDense,
2571: MatGetDiagonal_SeqDense,
2572: MatDiagonalScale_SeqDense,
2573: MatNorm_SeqDense,
2574: /* 20*/ MatAssemblyBegin_SeqDense,
2575: MatAssemblyEnd_SeqDense,
2576: MatSetOption_SeqDense,
2577: MatZeroEntries_SeqDense,
2578: /* 24*/ MatZeroRows_SeqDense,
2579: NULL,
2580: NULL,
2581: NULL,
2582: NULL,
2583: /* 29*/ MatSetUp_SeqDense,
2584: NULL,
2585: NULL,
2586: NULL,
2587: NULL,
2588: /* 34*/ MatDuplicate_SeqDense,
2589: NULL,
2590: NULL,
2591: NULL,
2592: NULL,
2593: /* 39*/ MatAXPY_SeqDense,
2594: MatCreateSubMatrices_SeqDense,
2595: NULL,
2596: MatGetValues_SeqDense,
2597: MatCopy_SeqDense,
2598: /* 44*/ MatGetRowMax_SeqDense,
2599: MatScale_SeqDense,
2600: MatShift_Basic,
2601: NULL,
2602: MatZeroRowsColumns_SeqDense,
2603: /* 49*/ MatSetRandom_SeqDense,
2604: NULL,
2605: NULL,
2606: NULL,
2607: NULL,
2608: /* 54*/ NULL,
2609: NULL,
2610: NULL,
2611: NULL,
2612: NULL,
2613: /* 59*/ NULL,
2614: MatDestroy_SeqDense,
2615: MatView_SeqDense,
2616: NULL,
2617: NULL,
2618: /* 64*/ NULL,
2619: NULL,
2620: NULL,
2621: NULL,
2622: NULL,
2623: /* 69*/ MatGetRowMaxAbs_SeqDense,
2624: NULL,
2625: NULL,
2626: NULL,
2627: NULL,
2628: /* 74*/ NULL,
2629: NULL,
2630: NULL,
2631: NULL,
2632: NULL,
2633: /* 79*/ NULL,
2634: NULL,
2635: NULL,
2636: NULL,
2637: /* 83*/ MatLoad_SeqDense,
2638: MatIsSymmetric_SeqDense,
2639: MatIsHermitian_SeqDense,
2640: NULL,
2641: NULL,
2642: NULL,
2643: /* 89*/ NULL,
2644: NULL,
2645: MatMatMultNumeric_SeqDense_SeqDense,
2646: NULL,
2647: NULL,
2648: /* 94*/ NULL,
2649: NULL,
2650: NULL,
2651: MatMatTransposeMultNumeric_SeqDense_SeqDense,
2652: NULL,
2653: /* 99*/ MatProductSetFromOptions_SeqDense,
2654: NULL,
2655: NULL,
2656: MatConjugate_SeqDense,
2657: NULL,
2658: /*104*/ NULL,
2659: MatRealPart_SeqDense,
2660: MatImaginaryPart_SeqDense,
2661: NULL,
2662: NULL,
2663: /*109*/ NULL,
2664: NULL,
2665: MatGetRowMin_SeqDense,
2666: MatGetColumnVector_SeqDense,
2667: MatMissingDiagonal_SeqDense,
2668: /*114*/ NULL,
2669: NULL,
2670: NULL,
2671: NULL,
2672: NULL,
2673: /*119*/ NULL,
2674: NULL,
2675: NULL,
2676: NULL,
2677: NULL,
2678: /*124*/ NULL,
2679: MatGetColumnNorms_SeqDense,
2680: NULL,
2681: NULL,
2682: NULL,
2683: /*129*/ NULL,
2684: NULL,
2685: NULL,
2686: MatTransposeMatMultNumeric_SeqDense_SeqDense,
2687: NULL,
2688: /*134*/ NULL,
2689: NULL,
2690: NULL,
2691: NULL,
2692: NULL,
2693: /*139*/ NULL,
2694: NULL,
2695: NULL,
2696: NULL,
2697: NULL,
2698: MatCreateMPIMatConcatenateSeqMat_SeqDense,
2699: /*145*/ NULL,
2700: NULL,
2701: NULL
2702: };
2704: /*@C
2705: MatCreateSeqDense - Creates a sequential dense matrix that
2706: is stored in column major order (the usual Fortran 77 manner). Many
2707: of the matrix operations use the BLAS and LAPACK routines.
2709: Collective
2711: Input Parameters:
2712: + comm - MPI communicator, set to PETSC_COMM_SELF
2713: . m - number of rows
2714: . n - number of columns
2715: - data - optional location of matrix data in column major order. Set data=NULL for PETSc
2716: to control all matrix memory allocation.
2718: Output Parameter:
2719: . A - the matrix
2721: Notes:
2722: The data input variable is intended primarily for Fortran programmers
2723: who wish to allocate their own matrix memory space. Most users should
2724: set data=NULL.
2726: Level: intermediate
2728: .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2729: @*/
2730: PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2731: {
2735: MatCreate(comm,A);
2736: MatSetSizes(*A,m,n,m,n);
2737: MatSetType(*A,MATSEQDENSE);
2738: MatSeqDenseSetPreallocation(*A,data);
2739: return(0);
2740: }
2742: /*@C
2743: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
2745: Collective
2747: Input Parameters:
2748: + B - the matrix
2749: - data - the array (or NULL)
2751: Notes:
2752: The data input variable is intended primarily for Fortran programmers
2753: who wish to allocate their own matrix memory space. Most users should
2754: need not call this routine.
2756: Level: intermediate
2758: .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatDenseSetLDA()
2760: @*/
2761: PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2762: {
2767: PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));
2768: return(0);
2769: }
2771: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2772: {
2773: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2777: if (b->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
2778: B->preallocated = PETSC_TRUE;
2780: PetscLayoutSetUp(B->rmap);
2781: PetscLayoutSetUp(B->cmap);
2783: if (b->lda <= 0) b->lda = B->rmap->n;
2785: PetscIntMultError(b->lda,B->cmap->n,NULL);
2786: if (!data) { /* petsc-allocated storage */
2787: if (!b->user_alloc) { PetscFree(b->v); }
2788: PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v);
2789: PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar));
2791: b->user_alloc = PETSC_FALSE;
2792: } else { /* user-allocated storage */
2793: if (!b->user_alloc) { PetscFree(b->v); }
2794: b->v = data;
2795: b->user_alloc = PETSC_TRUE;
2796: }
2797: B->assembled = PETSC_TRUE;
2798: return(0);
2799: }
2801: #if defined(PETSC_HAVE_ELEMENTAL)
2802: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
2803: {
2804: Mat mat_elemental;
2805: PetscErrorCode ierr;
2806: const PetscScalar *array;
2807: PetscScalar *v_colwise;
2808: PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols;
2811: PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);
2812: MatDenseGetArrayRead(A,&array);
2813: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
2814: k = 0;
2815: for (j=0; j<N; j++) {
2816: cols[j] = j;
2817: for (i=0; i<M; i++) {
2818: v_colwise[j*M+i] = array[k++];
2819: }
2820: }
2821: for (i=0; i<M; i++) {
2822: rows[i] = i;
2823: }
2824: MatDenseRestoreArrayRead(A,&array);
2826: MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);
2827: MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);
2828: MatSetType(mat_elemental,MATELEMENTAL);
2829: MatSetUp(mat_elemental);
2831: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2832: MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);
2833: MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);
2834: MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);
2835: PetscFree3(v_colwise,rows,cols);
2837: if (reuse == MAT_INPLACE_MATRIX) {
2838: MatHeaderReplace(A,&mat_elemental);
2839: } else {
2840: *newmat = mat_elemental;
2841: }
2842: return(0);
2843: }
2844: #endif
2846: PetscErrorCode MatDenseSetLDA_SeqDense(Mat B,PetscInt lda)
2847: {
2848: Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2849: PetscBool data;
2852: data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
2853: if (!b->user_alloc && data && b->lda!=lda) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"LDA cannot be changed after allocation of internal storage");
2854: 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);
2855: b->lda = lda;
2856: return(0);
2857: }
2859: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2860: {
2862: PetscMPIInt size;
2865: MPI_Comm_size(comm,&size);
2866: if (size == 1) {
2867: if (scall == MAT_INITIAL_MATRIX) {
2868: MatDuplicate(inmat,MAT_COPY_VALUES,outmat);
2869: } else {
2870: MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
2871: }
2872: } else {
2873: MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);
2874: }
2875: return(0);
2876: }
2878: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
2879: {
2880: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2884: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
2885: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
2886: if (!a->cvec) {
2887: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
2888: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
2889: }
2890: a->vecinuse = col + 1;
2891: MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);
2892: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
2893: *v = a->cvec;
2894: return(0);
2895: }
2897: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
2898: {
2899: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2903: if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
2904: if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
2905: a->vecinuse = 0;
2906: MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);
2907: VecResetArray(a->cvec);
2908: *v = NULL;
2909: return(0);
2910: }
2912: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
2913: {
2914: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2918: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
2919: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
2920: if (!a->cvec) {
2921: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
2922: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
2923: }
2924: a->vecinuse = col + 1;
2925: MatDenseGetArrayRead(A,&a->ptrinuse);
2926: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
2927: VecLockReadPush(a->cvec);
2928: *v = a->cvec;
2929: return(0);
2930: }
2932: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
2933: {
2934: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2938: if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
2939: if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
2940: a->vecinuse = 0;
2941: MatDenseRestoreArrayRead(A,&a->ptrinuse);
2942: VecLockReadPop(a->cvec);
2943: VecResetArray(a->cvec);
2944: *v = NULL;
2945: return(0);
2946: }
2948: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
2949: {
2950: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2954: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
2955: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
2956: if (!a->cvec) {
2957: VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);
2958: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);
2959: }
2960: a->vecinuse = col + 1;
2961: MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);
2962: VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);
2963: *v = a->cvec;
2964: return(0);
2965: }
2967: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
2968: {
2969: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2973: if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first");
2974: if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector");
2975: a->vecinuse = 0;
2976: MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);
2977: VecResetArray(a->cvec);
2978: *v = NULL;
2979: return(0);
2980: }
2982: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
2983: {
2984: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2988: if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first");
2989: if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
2990: if (a->cmat && cend-cbegin != a->cmat->cmap->N) {
2991: MatDestroy(&a->cmat);
2992: }
2993: if (!a->cmat) {
2994: 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);
2995: PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat);
2996: } else {
2997: MatDensePlaceArray(a->cmat,a->v + (size_t)cbegin * (size_t)a->lda);
2998: }
2999: MatDenseSetLDA(a->cmat,a->lda);
3000: a->matinuse = cbegin + 1;
3001: *v = a->cmat;
3002: return(0);
3003: }
3005: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v)
3006: {
3007: Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3011: if (!a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetSubMatrix() first");
3012: if (!a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column matrix");
3013: if (*v != a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not the matrix obtained from MatDenseGetSubMatrix()");
3014: a->matinuse = 0;
3015: MatDenseResetArray(a->cmat);
3016: *v = NULL;
3017: return(0);
3018: }
3020: /*MC
3021: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
3023: Options Database Keys:
3024: . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
3026: Level: beginner
3028: .seealso: MatCreateSeqDense()
3030: M*/
3031: PetscErrorCode MatCreate_SeqDense(Mat B)
3032: {
3033: Mat_SeqDense *b;
3035: PetscMPIInt size;
3038: MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
3039: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
3041: PetscNewLog(B,&b);
3042: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
3043: B->data = (void*)b;
3045: b->roworiented = PETSC_TRUE;
3047: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);
3048: PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense);
3049: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);
3050: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);
3051: PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);
3052: PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);
3053: PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense);
3054: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);
3055: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);
3056: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense);
3057: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense);
3058: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);
3059: #if defined(PETSC_HAVE_ELEMENTAL)
3060: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);
3061: #endif
3062: #if defined(PETSC_HAVE_SCALAPACK)
3063: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_scalapack_C",MatConvert_Dense_ScaLAPACK);
3064: #endif
3065: #if defined(PETSC_HAVE_CUDA)
3066: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);
3067: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3068: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);
3069: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);
3070: #endif
3071: PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);
3072: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);
3073: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);
3074: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
3075: PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);
3077: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);
3078: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);
3079: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);
3080: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);
3081: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);
3082: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);
3083: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);
3084: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);
3085: PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense);
3086: PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense);
3087: PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);
3088: return(0);
3089: }
3091: /*@C
3092: 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.
3094: Not Collective
3096: Input Parameters:
3097: + mat - a MATSEQDENSE or MATMPIDENSE matrix
3098: - col - column index
3100: Output Parameter:
3101: . vals - pointer to the data
3103: Level: intermediate
3105: .seealso: MatDenseRestoreColumn()
3106: @*/
3107: PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
3108: {
3115: PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));
3116: return(0);
3117: }
3119: /*@C
3120: MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn().
3122: Not Collective
3124: Input Parameter:
3125: . mat - a MATSEQDENSE or MATMPIDENSE matrix
3127: Output Parameter:
3128: . vals - pointer to the data
3130: Level: intermediate
3132: .seealso: MatDenseGetColumn()
3133: @*/
3134: PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
3135: {
3141: PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));
3142: return(0);
3143: }
3145: /*@C
3146: MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec.
3148: Collective
3150: Input Parameters:
3151: + mat - the Mat object
3152: - col - the column index
3154: Output Parameter:
3155: . v - the vector
3157: Notes:
3158: The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVec() when the vector is no longer needed.
3159: Use MatDenseGetColumnVecRead() to obtain read-only access or MatDenseGetColumnVecWrite() for write-only access.
3161: Level: intermediate
3163: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3164: @*/
3165: PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v)
3166: {
3174: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3175: 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);
3176: PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3177: return(0);
3178: }
3180: /*@C
3181: MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec().
3183: Collective
3185: Input Parameters:
3186: + mat - the Mat object
3187: . col - the column index
3188: - v - the Vec object
3190: Level: intermediate
3192: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3193: @*/
3194: PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v)
3195: {
3203: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3204: 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);
3205: PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));
3206: return(0);
3207: }
3209: /*@C
3210: MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec.
3212: Collective
3214: Input Parameters:
3215: + mat - the Mat object
3216: - col - the column index
3218: Output Parameter:
3219: . v - the vector
3221: Notes:
3222: The vector is owned by PETSc and users cannot modify it.
3223: Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed.
3224: Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access.
3226: Level: intermediate
3228: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3229: @*/
3230: PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v)
3231: {
3239: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3240: 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);
3241: PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3242: return(0);
3243: }
3245: /*@C
3246: MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead().
3248: Collective
3250: Input Parameters:
3251: + mat - the Mat object
3252: . col - the column index
3253: - v - the Vec object
3255: Level: intermediate
3257: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite()
3258: @*/
3259: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v)
3260: {
3268: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3269: 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);
3270: PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));
3271: return(0);
3272: }
3274: /*@C
3275: MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec.
3277: Collective
3279: Input Parameters:
3280: + mat - the Mat object
3281: - col - the column index
3283: Output Parameter:
3284: . v - the vector
3286: Notes:
3287: The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed.
3288: Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access.
3290: Level: intermediate
3292: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3293: @*/
3294: PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v)
3295: {
3303: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3304: 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);
3305: PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3306: return(0);
3307: }
3309: /*@C
3310: MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite().
3312: Collective
3314: Input Parameters:
3315: + mat - the Mat object
3316: . col - the column index
3317: - v - the Vec object
3319: Level: intermediate
3321: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead()
3322: @*/
3323: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v)
3324: {
3332: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3333: 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);
3334: PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));
3335: return(0);
3336: }
3338: /*@C
3339: MatDenseGetSubMatrix - Gives access to a block of columns of a dense matrix, represented as a Mat.
3341: Collective
3343: Input Parameters:
3344: + mat - the Mat object
3345: . cbegin - the first index in the block
3346: - cend - the last index in the block
3348: Output Parameter:
3349: . v - the matrix
3351: Notes:
3352: The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed.
3354: Level: intermediate
3356: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseRestoreSubMatrix()
3357: @*/
3358: PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3359: {
3368: if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated");
3369: 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);
3370: 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);
3371: PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));
3372: return(0);
3373: }
3375: /*@C
3376: MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix().
3378: Collective
3380: Input Parameters:
3381: + mat - the Mat object
3382: - v - the Mat object
3384: Level: intermediate
3386: .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseGetSubMatrix()
3387: @*/
3388: PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v)
3389: {
3396: PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v));
3397: return(0);
3398: }