1
2 /*
3 Defines the basic matrix operations for sequential dense.
4 */
5
6 #include <../src/mat/impls/dense/seq/dense.h> /*I "petscmat.h" I*/
7 #include <petscblaslapack.h>
8
9 #include <../src/mat/impls/aij/seq/aij.h>
10
MatSeqDenseSymmetrize_Private(Mat A,PetscBool hermitian)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;
16 PetscErrorCode ierr;
17
18 PetscFunctionBegin;
19 if (A->rmap->n != A->cmap->n) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot symmetrize a rectangular matrix");
20 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr);
35 PetscFunctionReturn(0);
36 }
37
MatSeqDenseInvertFactors_Private(Mat A)38 PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A)
39 {
40 Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
41 PetscErrorCode ierr;
42 PetscBLASInt info,n;
43
44 PetscFunctionBegin;
45 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0);
46 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr);
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 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr);
52 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr);
53 }
54 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
55 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
56 ierr = PetscFPTrapPop();CHKERRQ(ierr);
57 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr);
58 } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
59 if (A->spd) {
60 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
61 PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info));
62 ierr = PetscFPTrapPop();CHKERRQ(ierr);
63 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);CHKERRQ(ierr);
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 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
69 PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
70 ierr = PetscFPTrapPop();CHKERRQ(ierr);
71 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);CHKERRQ(ierr);
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 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
77 PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info));
78 ierr = PetscFPTrapPop();CHKERRQ(ierr);
79 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_FALSE);CHKERRQ(ierr);
80 }
81 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad Inversion: zero pivot in row %D",(PetscInt)info-1);
82 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr);
83 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve");
84
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 ierr = PetscFree(A->solvertype);CHKERRQ(ierr);
93 PetscFunctionReturn(0);
94 }
95
MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)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;
103
104 PetscFunctionBegin;
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) PetscFunctionReturn(0);
113
114 /* fix right hand side if needed */
115 if (x && b) {
116 Vec xt;
117
118 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices");
119 ierr = VecDuplicate(x,&xt);CHKERRQ(ierr);
120 ierr = VecCopy(x,xt);CHKERRQ(ierr);
121 ierr = VecScale(xt,-1.0);CHKERRQ(ierr);
122 ierr = MatMultAdd(A,xt,b,b);CHKERRQ(ierr);
123 ierr = VecDestroy(&xt);CHKERRQ(ierr);
124 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
125 ierr = VecGetArray(b,&bb);CHKERRQ(ierr);
126 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
127 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
128 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr);
129 }
130
131 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr);
132 for (i=0; i<N; i++) {
133 slot = v + rows[i]*m;
134 ierr = PetscArrayzero(slot,r);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr);
148 PetscFunctionReturn(0);
149 }
150
MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)151 PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C)
152 {
153 Mat_SeqDense *c = (Mat_SeqDense*)(C->data);
154 PetscErrorCode ierr;
155
156 PetscFunctionBegin;
157 if (c->ptapwork) {
158 ierr = (*C->ops->matmultnumeric)(A,P,c->ptapwork);CHKERRQ(ierr);
159 ierr = (*C->ops->transposematmultnumeric)(P,c->ptapwork,C);CHKERRQ(ierr);
160 } else SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Must call MatPtAPSymbolic_SeqDense_SeqDense() first");
161 PetscFunctionReturn(0);
162 }
163
MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat C)164 PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat C)
165 {
166 Mat_SeqDense *c;
167 PetscBool cisdense;
168 PetscErrorCode ierr;
169
170 PetscFunctionBegin;
171 ierr = MatSetSizes(C,P->cmap->n,P->cmap->n,P->cmap->N,P->cmap->N);CHKERRQ(ierr);
172 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr);
173 if (!cisdense) {
174 PetscBool flg;
175
176 ierr = PetscObjectTypeCompare((PetscObject)P,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr);
177 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr);
178 }
179 ierr = MatSetUp(C);CHKERRQ(ierr);
180 c = (Mat_SeqDense*)C->data;
181 ierr = MatCreate(PetscObjectComm((PetscObject)A),&c->ptapwork);CHKERRQ(ierr);
182 ierr = MatSetSizes(c->ptapwork,A->rmap->n,P->cmap->n,A->rmap->N,P->cmap->N);CHKERRQ(ierr);
183 ierr = MatSetType(c->ptapwork,((PetscObject)C)->type_name);CHKERRQ(ierr);
184 ierr = MatSetUp(c->ptapwork);CHKERRQ(ierr);
185 PetscFunctionReturn(0);
186 }
187
MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat * newmat)188 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
189 {
190 Mat B = NULL;
191 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
192 Mat_SeqDense *b;
193 PetscErrorCode ierr;
194 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i;
195 MatScalar *av=a->a;
196 PetscBool isseqdense;
197
198 PetscFunctionBegin;
199 if (reuse == MAT_REUSE_MATRIX) {
200 ierr = PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);CHKERRQ(ierr);
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 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr);
205 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr);
206 ierr = MatSetType(B,MATSEQDENSE);CHKERRQ(ierr);
207 ierr = MatSeqDenseSetPreallocation(B,NULL);CHKERRQ(ierr);
208 b = (Mat_SeqDense*)(B->data);
209 } else {
210 b = (Mat_SeqDense*)((*newmat)->data);
211 ierr = PetscArrayzero(b->v,m*n);CHKERRQ(ierr);
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 }
222
223 if (reuse == MAT_INPLACE_MATRIX) {
224 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
225 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
226 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr);
227 } else {
228 if (B) *newmat = B;
229 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
230 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
231 }
232 PetscFunctionReturn(0);
233 }
234
MatConvert_SeqDense_SeqAIJ(Mat A,MatType newtype,MatReuse reuse,Mat * newmat)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;
239 PetscErrorCode ierr;
240 PetscInt i, j;
241 PetscInt *rows, *nnz;
242 MatScalar *aa = a->v, *vals;
243
244 PetscFunctionBegin;
245 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr);
246 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr);
247 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr);
248 ierr = PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);CHKERRQ(ierr);
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 ierr = MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);CHKERRQ(ierr);
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 ierr = MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);CHKERRQ(ierr);
259 aa += a->lda;
260 }
261 ierr = PetscFree3(rows,nnz,vals);CHKERRQ(ierr);
262 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
263 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
264
265 if (reuse == MAT_INPLACE_MATRIX) {
266 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr);
267 } else {
268 *newmat = B;
269 }
270 PetscFunctionReturn(0);
271 }
272
MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str)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;
280
281 PetscFunctionBegin;
282 ierr = MatDenseGetArrayRead(X,&xv);CHKERRQ(ierr);
283 ierr = MatDenseGetArray(Y,&yv);CHKERRQ(ierr);
284 ierr = PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);CHKERRQ(ierr);
285 ierr = PetscBLASIntCast(X->rmap->n,&m);CHKERRQ(ierr);
286 ierr = PetscBLASIntCast(x->lda,&ldax);CHKERRQ(ierr);
287 ierr = PetscBLASIntCast(y->lda,&lday);CHKERRQ(ierr);
288 if (ldax>m || lday>m) {
289 PetscInt j;
290
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 ierr = MatDenseRestoreArrayRead(X,&xv);CHKERRQ(ierr);
298 ierr = MatDenseRestoreArray(Y,&yv);CHKERRQ(ierr);
299 ierr = PetscLogFlops(PetscMax(2.0*N-1,0));CHKERRQ(ierr);
300 PetscFunctionReturn(0);
301 }
302
MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo * info)303 static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info)
304 {
305 PetscLogDouble N = A->rmap->n*A->cmap->n;
306
307 PetscFunctionBegin;
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 PetscFunctionReturn(0);
319 }
320
MatScale_SeqDense(Mat A,PetscScalar alpha)321 PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha)
322 {
323 Mat_SeqDense *a = (Mat_SeqDense*)A->data;
324 PetscScalar *v;
325 PetscErrorCode ierr;
326 PetscBLASInt one = 1,j,nz,lda = 0;
327
328 PetscFunctionBegin;
329 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr);
330 ierr = PetscBLASIntCast(a->lda,&lda);CHKERRQ(ierr);
331 if (lda>A->rmap->n) {
332 ierr = PetscBLASIntCast(A->rmap->n,&nz);CHKERRQ(ierr);
333 for (j=0; j<A->cmap->n; j++) {
334 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v+j*lda,&one));
335 }
336 } else {
337 ierr = PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);CHKERRQ(ierr);
338 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v,&one));
339 }
340 ierr = PetscLogFlops(nz);CHKERRQ(ierr);
341 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr);
342 PetscFunctionReturn(0);
343 }
344
MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool * fl)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;
351
352 PetscFunctionBegin;
353 *fl = PETSC_FALSE;
354 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0);
355 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr);
366 PetscFunctionReturn(0);
367 }
368
MatIsSymmetric_SeqDense(Mat A,PetscReal rtol,PetscBool * fl)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;
375
376 PetscFunctionBegin;
377 *fl = PETSC_FALSE;
378 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0);
379 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr);
390 PetscFunctionReturn(0);
391 }
392
MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)393 PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues)
394 {
395 Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
396 PetscErrorCode ierr;
397 PetscInt lda = (PetscInt)mat->lda,j,m,nlda = lda;
398
399 PetscFunctionBegin;
400 ierr = PetscLayoutReference(A->rmap,&newi->rmap);CHKERRQ(ierr);
401 ierr = PetscLayoutReference(A->cmap,&newi->cmap);CHKERRQ(ierr);
402 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */
403 ierr = MatDenseSetLDA(newi,lda);CHKERRQ(ierr);
404 }
405 ierr = MatSeqDenseSetPreallocation(newi,NULL);CHKERRQ(ierr);
406 if (cpvalues == MAT_COPY_VALUES) {
407 const PetscScalar *av;
408 PetscScalar *v;
409
410 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr);
411 ierr = MatDenseGetArray(newi,&v);CHKERRQ(ierr);
412 ierr = MatDenseGetLDA(newi,&nlda);CHKERRQ(ierr);
413 m = A->rmap->n;
414 if (lda>m || nlda>m) {
415 for (j=0; j<A->cmap->n; j++) {
416 ierr = PetscArraycpy(v+j*nlda,av+j*lda,m);CHKERRQ(ierr);
417 }
418 } else {
419 ierr = PetscArraycpy(v,av,A->rmap->n*A->cmap->n);CHKERRQ(ierr);
420 }
421 ierr = MatDenseRestoreArray(newi,&v);CHKERRQ(ierr);
422 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr);
423 }
424 PetscFunctionReturn(0);
425 }
426
MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat * newmat)427 PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat)
428 {
429 PetscErrorCode ierr;
430
431 PetscFunctionBegin;
432 ierr = MatCreate(PetscObjectComm((PetscObject)A),newmat);CHKERRQ(ierr);
433 ierr = MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr);
434 ierr = MatSetType(*newmat,((PetscObject)A)->type_name);CHKERRQ(ierr);
435 ierr = MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);CHKERRQ(ierr);
436 PetscFunctionReturn(0);
437 }
438
MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo * info_dummy)439 static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
440 {
441 MatFactorInfo info;
442 PetscErrorCode ierr;
443
444 PetscFunctionBegin;
445 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr);
446 ierr = (*fact->ops->lufactor)(fact,NULL,NULL,&info);CHKERRQ(ierr);
447 PetscFunctionReturn(0);
448 }
449
MatSolve_SeqDense(Mat A,Vec xx,Vec yy)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;
457
458 PetscFunctionBegin;
459 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr);
460 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr);
461 ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
462 ierr = PetscArraycpy(y,x,A->rmap->n);CHKERRQ(ierr);
463 if (A->factortype == MAT_FACTOR_LU) {
464 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
465 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
466 ierr = PetscFPTrapPop();CHKERRQ(ierr);
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 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
471 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
472 ierr = PetscFPTrapPop();CHKERRQ(ierr);
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 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
477 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
478 ierr = PetscFPTrapPop();CHKERRQ(ierr);
479 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
480 #endif
481 } else { /* symmetric case */
482 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
483 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
484 ierr = PetscFPTrapPop();CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr);
489 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr);
490 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr);
491 PetscFunctionReturn(0);
492 }
493
MatMatSolve_SeqDense(Mat A,Mat B,Mat X)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;
502
503 PetscFunctionBegin;
504 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr);
505 ierr = MatGetSize(B,NULL,&n);CHKERRQ(ierr);
506 ierr = PetscBLASIntCast(n,&nrhs);CHKERRQ(ierr);
507 ierr = MatDenseGetArrayRead(B,&b);CHKERRQ(ierr);
508 ierr = MatDenseGetArray(X,&x);CHKERRQ(ierr);
509
510 ierr = PetscArraycpy(x,b,m*nrhs);CHKERRQ(ierr);
511
512 if (A->factortype == MAT_FACTOR_LU) {
513 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
514 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
515 ierr = PetscFPTrapPop();CHKERRQ(ierr);
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 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
520 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info));
521 ierr = PetscFPTrapPop();CHKERRQ(ierr);
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 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
526 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
527 ierr = PetscFPTrapPop();CHKERRQ(ierr);
528 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve");
529 #endif
530 } else { /* symmetric case */
531 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
532 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info));
533 ierr = PetscFPTrapPop();CHKERRQ(ierr);
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");
537
538 ierr = MatDenseRestoreArrayRead(B,&b);CHKERRQ(ierr);
539 ierr = MatDenseRestoreArray(X,&x);CHKERRQ(ierr);
540 ierr = PetscLogFlops(nrhs*(2.0*m*m - m));CHKERRQ(ierr);
541 PetscFunctionReturn(0);
542 }
543
544 static PetscErrorCode MatConjugate_SeqDense(Mat);
545
MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy)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;
553
554 PetscFunctionBegin;
555 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr);
556 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr);
557 ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
558 ierr = PetscArraycpy(y,x,A->rmap->n);CHKERRQ(ierr);
559 if (A->factortype == MAT_FACTOR_LU) {
560 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
561 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
562 ierr = PetscFPTrapPop();CHKERRQ(ierr);
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 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr);
568 #endif
569 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
570 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info));
571 ierr = PetscFPTrapPop();CHKERRQ(ierr);
572 #if defined(PETSC_USE_COMPLEX)
573 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr);
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 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr);
579 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
580 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
581 ierr = PetscFPTrapPop();CHKERRQ(ierr);
582 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr);
583 #endif
584 } else { /* symmetric case */
585 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
586 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info));
587 ierr = PetscFPTrapPop();CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr);
592 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr);
593 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr);
594 PetscFunctionReturn(0);
595 }
596
597 /* ---------------------------------------------------------------*/
598 /* COMMENT: I have chosen to hide row permutation in the pivots,
599 rather than put it in the Mat->row slot.*/
MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo * minfo)600 PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo)
601 {
602 Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
603 PetscErrorCode ierr;
604 PetscBLASInt n,m,info;
605
606 PetscFunctionBegin;
607 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr);
608 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr);
609 if (!mat->pivots) {
610 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr);
611 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr);
612 }
613 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0);
614 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
615 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info));
616 ierr = PetscFPTrapPop();CHKERRQ(ierr);
617
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");
620
621 A->ops->solve = MatSolve_SeqDense;
622 A->ops->matsolve = MatMatSolve_SeqDense;
623 A->ops->solvetranspose = MatSolveTranspose_SeqDense;
624 A->factortype = MAT_FACTOR_LU;
625
626 ierr = PetscFree(A->solvertype);CHKERRQ(ierr);
627 ierr = PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);CHKERRQ(ierr);
628
629 ierr = PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);CHKERRQ(ierr);
630 PetscFunctionReturn(0);
631 }
632
633 /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo * factinfo)634 PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo)
635 {
636 Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
637 PetscErrorCode ierr;
638 PetscBLASInt info,n;
639
640 PetscFunctionBegin;
641 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr);
642 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0);
643 if (A->spd) {
644 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
645 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info));
646 ierr = PetscFPTrapPop();CHKERRQ(ierr);
647 #if defined(PETSC_USE_COMPLEX)
648 } else if (A->hermitian) {
649 if (!mat->pivots) {
650 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr);
651 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr);
652 }
653 if (!mat->fwork) {
654 PetscScalar dummy;
655
656 mat->lfwork = -1;
657 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
658 PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
659 ierr = PetscFPTrapPop();CHKERRQ(ierr);
660 mat->lfwork = (PetscInt)PetscRealPart(dummy);
661 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr);
662 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr);
663 }
664 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
665 PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
666 ierr = PetscFPTrapPop();CHKERRQ(ierr);
667 #endif
668 } else { /* symmetric case */
669 if (!mat->pivots) {
670 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr);
671 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr);
672 }
673 if (!mat->fwork) {
674 PetscScalar dummy;
675
676 mat->lfwork = -1;
677 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
678 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info));
679 ierr = PetscFPTrapPop();CHKERRQ(ierr);
680 mat->lfwork = (PetscInt)PetscRealPart(dummy);
681 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr);
682 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr);
683 }
684 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
685 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info));
686 ierr = PetscFPTrapPop();CHKERRQ(ierr);
687 }
688 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1);
689
690 A->ops->solve = MatSolve_SeqDense;
691 A->ops->matsolve = MatMatSolve_SeqDense;
692 A->ops->solvetranspose = MatSolveTranspose_SeqDense;
693 A->factortype = MAT_FACTOR_CHOLESKY;
694
695 ierr = PetscFree(A->solvertype);CHKERRQ(ierr);
696 ierr = PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);CHKERRQ(ierr);
697
698 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr);
699 PetscFunctionReturn(0);
700 }
701
MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo * info_dummy)702 PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy)
703 {
704 PetscErrorCode ierr;
705 MatFactorInfo info;
706
707 PetscFunctionBegin;
708 info.fill = 1.0;
709
710 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr);
711 ierr = (*fact->ops->choleskyfactor)(fact,NULL,&info);CHKERRQ(ierr);
712 PetscFunctionReturn(0);
713 }
714
MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo * info)715 PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info)
716 {
717 PetscFunctionBegin;
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 PetscFunctionReturn(0);
725 }
726
MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo * info)727 PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info)
728 {
729 PetscFunctionBegin;
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 PetscFunctionReturn(0);
737 }
738
739 /* uses LAPACK */
MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat * fact)740 PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact)
741 {
742 PetscErrorCode ierr;
743
744 PetscFunctionBegin;
745 ierr = MatCreate(PetscObjectComm((PetscObject)A),fact);CHKERRQ(ierr);
746 ierr = MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr);
747 ierr = MatSetType(*fact,MATDENSE);CHKERRQ(ierr);
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;
755
756 ierr = PetscFree((*fact)->solvertype);CHKERRQ(ierr);
757 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);CHKERRQ(ierr);
758 PetscFunctionReturn(0);
759 }
760
761 /* ------------------------------------------------------------------*/
MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx)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;
770
771 PetscFunctionBegin;
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 ierr = PetscBLASIntCast(m,&bm);CHKERRQ(ierr);
777 if (flag & SOR_ZERO_INITIAL_GUESS) {
778 /* this is a hack fix, should have another version without the second BLASdotu */
779 ierr = VecSet(xx,zero);CHKERRQ(ierr);
780 }
781 ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
782 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr);
800 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
801 PetscFunctionReturn(0);
802 }
803
804 /* -----------------------------------------------------------------*/
MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy)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;
813
814 PetscFunctionBegin;
815 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr);
816 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr);
817 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr);
818 ierr = VecGetArrayWrite(yy,&y);CHKERRQ(ierr);
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 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr);
825 }
826 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr);
827 ierr = VecRestoreArrayWrite(yy,&y);CHKERRQ(ierr);
828 PetscFunctionReturn(0);
829 }
830
MatMult_SeqDense(Mat A,Vec xx,Vec yy)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;
838
839 PetscFunctionBegin;
840 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr);
841 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr);
842 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr);
843 ierr = VecGetArrayWrite(yy,&y);CHKERRQ(ierr);
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 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr);
850 }
851 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr);
852 ierr = VecRestoreArrayWrite(yy,&y);CHKERRQ(ierr);
853 PetscFunctionReturn(0);
854 }
855
MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)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;
863
864 PetscFunctionBegin;
865 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr);
866 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr);
867 ierr = VecCopy(zz,yy);CHKERRQ(ierr);
868 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0);
869 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr);
870 ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
871 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
872 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr);
873 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr);
874 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr);
875 PetscFunctionReturn(0);
876 }
877
MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy)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;
886
887 PetscFunctionBegin;
888 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr);
889 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr);
890 ierr = VecCopy(zz,yy);CHKERRQ(ierr);
891 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0);
892 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr);
893 ierr = VecGetArray(yy,&y);CHKERRQ(ierr);
894 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One));
895 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr);
896 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr);
897 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr);
898 PetscFunctionReturn(0);
899 }
900
901 /* -----------------------------------------------------------------*/
MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt * ncols,PetscInt ** cols,PetscScalar ** vals)902 static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
903 {
904 Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
905 PetscErrorCode ierr;
906 PetscInt i;
907
908 PetscFunctionBegin;
909 *ncols = A->cmap->n;
910 if (cols) {
911 ierr = PetscMalloc1(A->cmap->n+1,cols);CHKERRQ(ierr);
912 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i;
913 }
914 if (vals) {
915 const PetscScalar *v;
916
917 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr);
918 ierr = PetscMalloc1(A->cmap->n+1,vals);CHKERRQ(ierr);
919 v += row;
920 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;}
921 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr);
922 }
923 PetscFunctionReturn(0);
924 }
925
MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt * ncols,PetscInt ** cols,PetscScalar ** vals)926 static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
927 {
928 PetscErrorCode ierr;
929
930 PetscFunctionBegin;
931 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);}
932 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); }
933 PetscFunctionReturn(0);
934 }
935 /* ----------------------------------------------------------------*/
MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv)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;
945
946 PetscFunctionBegin;
947 ierr = MatDenseGetArray(A,&av);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArray(A,&av);CHKERRQ(ierr);
999 #if defined(PETSC_HAVE_CUDA)
1000 A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1001 #endif
1002 PetscFunctionReturn(0);
1003 }
1004
MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[])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;
1011
1012 PetscFunctionBegin;
1013 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr);
1025 PetscFunctionReturn(0);
1026 }
1027
1028 /* -----------------------------------------------------------------*/
1029
MatView_Dense_Binary(Mat mat,PetscViewer viewer)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;
1038
1039 PetscFunctionBegin;
1040 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr);
1041 ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);CHKERRQ(ierr);
1042 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
1043 if (skipHeader) format = PETSC_VIEWER_NATIVE;
1044
1045 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr);
1046
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) {ierr = PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);CHKERRQ(ierr);}
1051
1052 ierr = MatGetLocalSize(mat,&m,NULL);CHKERRQ(ierr);
1053 if (format != PETSC_VIEWER_NATIVE) {
1054 PetscInt nnz = m*N, *iwork;
1055 /* store row lengths for each row */
1056 ierr = PetscMalloc1(nnz,&iwork);CHKERRQ(ierr);
1057 for (i=0; i<m; i++) iwork[i] = N;
1058 ierr = PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr);
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 ierr = PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr);
1064 ierr = PetscFree(iwork);CHKERRQ(ierr);
1065 }
1066 /* store matrix values as a dense matrix in row major order */
1067 ierr = PetscMalloc1(m*N,&vwork);CHKERRQ(ierr);
1068 ierr = MatDenseGetArrayRead(mat,&v);CHKERRQ(ierr);
1069 ierr = MatDenseGetLDA(mat,&lda);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArrayRead(mat,&v);CHKERRQ(ierr);
1074 ierr = PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr);
1075 ierr = PetscFree(vwork);CHKERRQ(ierr);
1076 PetscFunctionReturn(0);
1077 }
1078
MatLoad_Dense_Binary(Mat mat,PetscViewer viewer)1079 PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer)
1080 {
1081 PetscErrorCode ierr;
1082 PetscBool skipHeader;
1083 PetscInt header[4],M,N,m,nz,lda,i,j,k;
1084 PetscInt rows,cols;
1085 PetscScalar *v,*vwork;
1086
1087 PetscFunctionBegin;
1088 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr);
1089 ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);CHKERRQ(ierr);
1090
1091 if (!skipHeader) {
1092 ierr = PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);CHKERRQ(ierr);
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 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr);
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 }
1104
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 ierr = MatSetUp(mat);CHKERRQ(ierr);
1109 /* check if global sizes are correct */
1110 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr);
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);
1112
1113 ierr = MatGetSize(mat,NULL,&N);CHKERRQ(ierr);
1114 ierr = MatGetLocalSize(mat,&m,NULL);CHKERRQ(ierr);
1115 ierr = MatDenseGetArray(mat,&v);CHKERRQ(ierr);
1116 ierr = MatDenseGetLDA(mat,&lda);CHKERRQ(ierr);
1117 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */
1118 PetscInt nnz = m*N;
1119 /* read in matrix values */
1120 ierr = PetscMalloc1(nnz,&vwork);CHKERRQ(ierr);
1121 ierr = PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr);
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 ierr = PetscFree(vwork);CHKERRQ(ierr);
1127 } else { /* matrix in file is sparse format */
1128 PetscInt nnz = 0, *rlens, *icols;
1129 /* read in row lengths */
1130 ierr = PetscMalloc1(m,&rlens);CHKERRQ(ierr);
1131 ierr = PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr);
1132 for (i=0; i<m; i++) nnz += rlens[i];
1133 /* read in column indices and values */
1134 ierr = PetscMalloc2(nnz,&icols,nnz,&vwork);CHKERRQ(ierr);
1135 ierr = PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr);
1136 ierr = PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr);
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 ierr = PetscFree(rlens);CHKERRQ(ierr);
1142 ierr = PetscFree2(icols,vwork);CHKERRQ(ierr);
1143 }
1144 ierr = MatDenseRestoreArray(mat,&v);CHKERRQ(ierr);
1145 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1146 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1147 PetscFunctionReturn(0);
1148 }
1149
MatLoad_SeqDense(Mat newMat,PetscViewer viewer)1150 PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1151 {
1152 PetscBool isbinary, ishdf5;
1153 PetscErrorCode ierr;
1154
1155 PetscFunctionBegin;
1156 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1);
1157 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
1158 /* force binary viewer to load .info file if it has not yet done so */
1159 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr);
1160 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
1161 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr);
1162 if (isbinary) {
1163 ierr = MatLoad_Dense_Binary(newMat,viewer);CHKERRQ(ierr);
1164 } else if (ishdf5) {
1165 #if defined(PETSC_HAVE_HDF5)
1166 ierr = MatLoad_Dense_HDF5(newMat,viewer);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
1174 }
1175
MatView_SeqDense_ASCII(Mat A,PetscViewer viewer)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
1187
1188 PetscFunctionBegin;
1189 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&av);CHKERRQ(ierr);
1190 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
1191 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1192 PetscFunctionReturn(0); /* do nothing for now */
1193 } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1194 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);
1195 for (i=0; i<A->rmap->n; i++) {
1196 v = av + i;
1197 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr);
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 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr);
1202 } else if (PetscRealPart(*v)) {
1203 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));CHKERRQ(ierr);
1204 }
1205 #else
1206 if (*v) {
1207 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);CHKERRQ(ierr);
1208 }
1209 #endif
1210 v += a->lda;
1211 }
1212 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
1213 }
1214 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);
1215 } else {
1216 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);
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 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr);
1226 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr);
1227 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr);
1228 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr);
1229 }
1230
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 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));CHKERRQ(ierr);
1237 } else {
1238 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr);
1239 }
1240 #else
1241 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);CHKERRQ(ierr);
1242 #endif
1243 v += a->lda;
1244 }
1245 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
1246 }
1247 if (format == PETSC_VIEWER_ASCII_MATLAB) {
1248 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr);
1249 }
1250 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);
1251 }
1252 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);CHKERRQ(ierr);
1253 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
1254 PetscFunctionReturn(0);
1255 }
1256
1257 #include <petscdraw.h>
MatView_SeqDense_Draw_Zoom(PetscDraw draw,void * Aa)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;
1268
1269 PetscFunctionBegin;
1270 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr);
1271 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
1272 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr);
1273
1274 /* Loop over matrix elements drawing boxes */
1275 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr);
1276 if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1277 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr);
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 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr);
1288 }
1289 }
1290 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr);
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;
1296
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 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr);
1302 ierr = PetscDrawScalePopup(popup,minv,maxv);CHKERRQ(ierr);
1303
1304 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr);
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 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr);
1313 }
1314 }
1315 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr);
1316 }
1317 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr);
1318 PetscFunctionReturn(0);
1319 }
1320
MatView_SeqDense_Draw(Mat A,PetscViewer viewer)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;
1326 PetscErrorCode ierr;
1327
1328 PetscFunctionBegin;
1329 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
1330 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr);
1331 if (isnull) PetscFunctionReturn(0);
1332
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 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr);
1336 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr);
1337 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr);
1338 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr);
1339 ierr = PetscDrawSave(draw);CHKERRQ(ierr);
1340 PetscFunctionReturn(0);
1341 }
1342
MatView_SeqDense(Mat A,PetscViewer viewer)1343 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer)
1344 {
1345 PetscErrorCode ierr;
1346 PetscBool iascii,isbinary,isdraw;
1347
1348 PetscFunctionBegin;
1349 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
1350 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
1351 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
1352
1353 if (iascii) {
1354 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr);
1355 } else if (isbinary) {
1356 ierr = MatView_Dense_Binary(A,viewer);CHKERRQ(ierr);
1357 } else if (isdraw) {
1358 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr);
1359 }
1360 PetscFunctionReturn(0);
1361 }
1362
MatDensePlaceArray_SeqDense(Mat A,const PetscScalar * array)1363 static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar *array)
1364 {
1365 Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1366
1367 PetscFunctionBegin;
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 PetscFunctionReturn(0);
1379 }
1380
MatDenseResetArray_SeqDense(Mat A)1381 static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1382 {
1383 Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1384
1385 PetscFunctionBegin;
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 PetscFunctionReturn(0);
1395 }
1396
MatDenseReplaceArray_SeqDense(Mat A,const PetscScalar * array)1397 static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A,const PetscScalar *array)
1398 {
1399 Mat_SeqDense *a = (Mat_SeqDense*)A->data;
1400 PetscErrorCode ierr;
1401
1402 PetscFunctionBegin;
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) { ierr = PetscFree(a->v);CHKERRQ(ierr); }
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 PetscFunctionReturn(0);
1412 }
1413
MatDestroy_SeqDense(Mat mat)1414 PetscErrorCode MatDestroy_SeqDense(Mat mat)
1415 {
1416 Mat_SeqDense *l = (Mat_SeqDense*)mat->data;
1417 PetscErrorCode ierr;
1418
1419 PetscFunctionBegin;
1420 #if defined(PETSC_USE_LOG)
1421 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n);
1422 #endif
1423 ierr = PetscFree(l->pivots);CHKERRQ(ierr);
1424 ierr = PetscFree(l->fwork);CHKERRQ(ierr);
1425 ierr = MatDestroy(&l->ptapwork);CHKERRQ(ierr);
1426 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);}
1427 if (!l->unplaced_user_alloc) {ierr = PetscFree(l->unplacedarray);CHKERRQ(ierr);}
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 ierr = VecDestroy(&l->cvec);CHKERRQ(ierr);
1431 ierr = MatDestroy(&l->cmat);CHKERRQ(ierr);
1432 ierr = PetscFree(mat->data);CHKERRQ(ierr);
1433
1434 ierr = PetscObjectChangeTypeName((PetscObject)mat,NULL);CHKERRQ(ierr);
1435 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);CHKERRQ(ierr);
1436 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseSetLDA_C",NULL);CHKERRQ(ierr);
1437 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);CHKERRQ(ierr);
1438 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr);
1439 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);CHKERRQ(ierr);
1440 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);CHKERRQ(ierr);
1441 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseReplaceArray_C",NULL);CHKERRQ(ierr);
1442 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);CHKERRQ(ierr);
1443 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);CHKERRQ(ierr);
1444 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayWrite_C",NULL);CHKERRQ(ierr);
1445 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayWrite_C",NULL);CHKERRQ(ierr);
1446 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);CHKERRQ(ierr);
1447 #if defined(PETSC_HAVE_ELEMENTAL)
1448 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);CHKERRQ(ierr);
1449 #endif
1450 #if defined(PETSC_HAVE_SCALAPACK)
1451 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_scalapack_C",NULL);CHKERRQ(ierr);
1452 #endif
1453 #if defined(PETSC_HAVE_CUDA)
1454 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);CHKERRQ(ierr);
1455 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL);CHKERRQ(ierr);
1456 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL);CHKERRQ(ierr);
1457 #endif
1458 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr);
1459 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL);CHKERRQ(ierr);
1460 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL);CHKERRQ(ierr);
1461 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL);CHKERRQ(ierr);
1462 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL);CHKERRQ(ierr);
1463
1464 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);CHKERRQ(ierr);
1465 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);CHKERRQ(ierr);
1466 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVec_C",NULL);CHKERRQ(ierr);
1467 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVec_C",NULL);CHKERRQ(ierr);
1468 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecRead_C",NULL);CHKERRQ(ierr);
1469 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecRead_C",NULL);CHKERRQ(ierr);
1470 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecWrite_C",NULL);CHKERRQ(ierr);
1471 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecWrite_C",NULL);CHKERRQ(ierr);
1472 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetSubMatrix_C",NULL);CHKERRQ(ierr);
1473 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreSubMatrix_C",NULL);CHKERRQ(ierr);
1474 PetscFunctionReturn(0);
1475 }
1476
MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat * matout)1477 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout)
1478 {
1479 Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1480 PetscErrorCode ierr;
1481 PetscInt k,j,m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1482 PetscScalar *v,tmp;
1483
1484 PetscFunctionBegin;
1485 if (reuse == MAT_INPLACE_MATRIX) {
1486 if (m == n) { /* in place transpose */
1487 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr);
1496 } else { /* reuse memory, temporary allocates new memory */
1497 PetscScalar *v2;
1498 PetscLayout tmplayout;
1499
1500 ierr = PetscMalloc1((size_t)m*n,&v2);CHKERRQ(ierr);
1501 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr);
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 ierr = PetscArraycpy(v,v2,(size_t)m*n);CHKERRQ(ierr);
1506 ierr = PetscFree(v2);CHKERRQ(ierr);
1507 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr);
1508 /* cleanup size dependent quantities */
1509 ierr = VecDestroy(&mat->cvec);CHKERRQ(ierr);
1510 ierr = MatDestroy(&mat->cmat);CHKERRQ(ierr);
1511 ierr = PetscFree(mat->pivots);CHKERRQ(ierr);
1512 ierr = PetscFree(mat->fwork);CHKERRQ(ierr);
1513 ierr = MatDestroy(&mat->ptapwork);CHKERRQ(ierr);
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;
1525
1526 if (reuse == MAT_INITIAL_MATRIX) {
1527 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr);
1528 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr);
1529 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr);
1530 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr);
1531 } else tmat = *matout;
1532
1533 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr);
1534 ierr = MatDenseGetArray(tmat,&v2);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArray(tmat,&v2);CHKERRQ(ierr);
1541 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr);
1542 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1543 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1544 *matout = tmat;
1545 }
1546 PetscFunctionReturn(0);
1547 }
1548
MatEqual_SeqDense(Mat A1,Mat A2,PetscBool * flg)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;
1556
1557 PetscFunctionBegin;
1558 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);}
1559 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);}
1560 ierr = MatDenseGetArrayRead(A1,&v1);CHKERRQ(ierr);
1561 ierr = MatDenseGetArrayRead(A2,&v2);CHKERRQ(ierr);
1562 for (i=0; i<A1->cmap->n; i++) {
1563 ierr = PetscArraycmp(v1,v2,A1->rmap->n,flg);CHKERRQ(ierr);
1564 if (*flg == PETSC_FALSE) PetscFunctionReturn(0);
1565 v1 += mat1->lda;
1566 v2 += mat2->lda;
1567 }
1568 ierr = MatDenseRestoreArrayRead(A1,&v1);CHKERRQ(ierr);
1569 ierr = MatDenseRestoreArrayRead(A2,&v2);CHKERRQ(ierr);
1570 *flg = PETSC_TRUE;
1571 PetscFunctionReturn(0);
1572 }
1573
MatGetDiagonal_SeqDense(Mat A,Vec v)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;
1581
1582 PetscFunctionBegin;
1583 ierr = VecGetSize(v,&n);CHKERRQ(ierr);
1584 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
1585 len = PetscMin(A->rmap->n,A->cmap->n);
1586 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr);
1592 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
1593 PetscFunctionReturn(0);
1594 }
1595
MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr)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;
1603
1604 PetscFunctionBegin;
1605 ierr = MatDenseGetArray(A,&vv);CHKERRQ(ierr);
1606 if (ll) {
1607 ierr = VecGetSize(ll,&m);CHKERRQ(ierr);
1608 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr);
1616 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr);
1617 }
1618 if (rr) {
1619 ierr = VecGetSize(rr,&n);CHKERRQ(ierr);
1620 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr);
1628 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr);
1629 }
1630 ierr = MatDenseRestoreArray(A,&vv);CHKERRQ(ierr);
1631 PetscFunctionReturn(0);
1632 }
1633
MatNorm_SeqDense(Mat A,NormType type,PetscReal * nrm)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;
1641
1642 PetscFunctionBegin;
1643 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr);
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 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr);
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 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr);
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 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr);
1688 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
1689 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr);
1690 PetscFunctionReturn(0);
1691 }
1692
MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)1693 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg)
1694 {
1695 Mat_SeqDense *aij = (Mat_SeqDense*)A->data;
1696 PetscErrorCode ierr;
1697
1698 PetscFunctionBegin;
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 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
1726 }
1727
MatZeroEntries_SeqDense(Mat A)1728 static PetscErrorCode MatZeroEntries_SeqDense(Mat A)
1729 {
1730 Mat_SeqDense *l = (Mat_SeqDense*)A->data;
1731 PetscErrorCode ierr;
1732 PetscInt lda=l->lda,m=A->rmap->n,j;
1733 PetscScalar *v;
1734
1735 PetscFunctionBegin;
1736 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr);
1737 if (lda>m) {
1738 for (j=0; j<A->cmap->n; j++) {
1739 ierr = PetscArrayzero(v+j*lda,m);CHKERRQ(ierr);
1740 }
1741 } else {
1742 ierr = PetscArrayzero(v,A->rmap->n*A->cmap->n);CHKERRQ(ierr);
1743 }
1744 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr);
1745 PetscFunctionReturn(0);
1746 }
1747
MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)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;
1755
1756 PetscFunctionBegin;
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) PetscFunctionReturn(0);
1764
1765 /* fix right hand side if needed */
1766 if (x && b) {
1767 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
1768 ierr = VecGetArray(b,&bb);CHKERRQ(ierr);
1769 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]];
1770 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
1771 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr);
1772 }
1773
1774 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr);
1787 PetscFunctionReturn(0);
1788 }
1789
MatDenseGetLDA_SeqDense(Mat A,PetscInt * lda)1790 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda)
1791 {
1792 Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1793
1794 PetscFunctionBegin;
1795 *lda = mat->lda;
1796 PetscFunctionReturn(0);
1797 }
1798
MatDenseGetArray_SeqDense(Mat A,PetscScalar ** array)1799 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array)
1800 {
1801 Mat_SeqDense *mat = (Mat_SeqDense*)A->data;
1802
1803 PetscFunctionBegin;
1804 if (mat->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
1805 *array = mat->v;
1806 PetscFunctionReturn(0);
1807 }
1808
MatDenseRestoreArray_SeqDense(Mat A,PetscScalar ** array)1809 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array)
1810 {
1811 PetscFunctionBegin;
1812 *array = NULL;
1813 PetscFunctionReturn(0);
1814 }
1815
1816 /*@C
1817 MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray()
1818
1819 Not collective
1820
1821 Input Parameter:
1822 . mat - a MATSEQDENSE or MATMPIDENSE matrix
1823
1824 Output Parameter:
1825 . lda - the leading dimension
1826
1827 Level: intermediate
1828
1829 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseSetLDA()
1830 @*/
MatDenseGetLDA(Mat A,PetscInt * lda)1831 PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda)
1832 {
1833 PetscErrorCode ierr;
1834
1835 PetscFunctionBegin;
1836 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
1837 PetscValidPointer(lda,2);
1838 ierr = PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));CHKERRQ(ierr);
1839 PetscFunctionReturn(0);
1840 }
1841
1842 /*@C
1843 MatDenseSetLDA - Sets the leading dimension of the array used by the dense matrix
1844
1845 Not collective
1846
1847 Input Parameter:
1848 + mat - a MATSEQDENSE or MATMPIDENSE matrix
1849 - lda - the leading dimension
1850
1851 Level: intermediate
1852
1853 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetLDA()
1854 @*/
MatDenseSetLDA(Mat A,PetscInt lda)1855 PetscErrorCode MatDenseSetLDA(Mat A,PetscInt lda)
1856 {
1857 PetscErrorCode ierr;
1858
1859 PetscFunctionBegin;
1860 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
1861 ierr = PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda));CHKERRQ(ierr);
1862 PetscFunctionReturn(0);
1863 }
1864
1865 /*@C
1866 MatDenseGetArray - gives read-write access to the array where the data for a dense matrix is stored
1867
1868 Logically Collective on Mat
1869
1870 Input Parameter:
1871 . mat - a dense matrix
1872
1873 Output Parameter:
1874 . array - pointer to the data
1875
1876 Level: intermediate
1877
1878 .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
1879 @*/
MatDenseGetArray(Mat A,PetscScalar ** array)1880 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array)
1881 {
1882 PetscErrorCode ierr;
1883
1884 PetscFunctionBegin;
1885 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
1886 PetscValidPointer(array,2);
1887 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr);
1888 PetscFunctionReturn(0);
1889 }
1890
1891 /*@C
1892 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray()
1893
1894 Logically Collective on Mat
1895
1896 Input Parameters:
1897 + mat - a dense matrix
1898 - array - pointer to the data
1899
1900 Level: intermediate
1901
1902 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
1903 @*/
MatDenseRestoreArray(Mat A,PetscScalar ** array)1904 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array)
1905 {
1906 PetscErrorCode ierr;
1907
1908 PetscFunctionBegin;
1909 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
1910 PetscValidPointer(array,2);
1911 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr);
1912 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr);
1913 #if defined(PETSC_HAVE_CUDA)
1914 A->offloadmask = PETSC_OFFLOAD_CPU;
1915 #endif
1916 PetscFunctionReturn(0);
1917 }
1918
1919 /*@C
1920 MatDenseGetArrayRead - gives read-only access to the array where the data for a dense matrix is stored
1921
1922 Not Collective
1923
1924 Input Parameter:
1925 . mat - a dense matrix
1926
1927 Output Parameter:
1928 . array - pointer to the data
1929
1930 Level: intermediate
1931
1932 .seealso: MatDenseRestoreArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
1933 @*/
MatDenseGetArrayRead(Mat A,const PetscScalar ** array)1934 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array)
1935 {
1936 PetscErrorCode ierr;
1937
1938 PetscFunctionBegin;
1939 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
1940 PetscValidPointer(array,2);
1941 ierr = PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr);
1942 PetscFunctionReturn(0);
1943 }
1944
1945 /*@C
1946 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayRead()
1947
1948 Not Collective
1949
1950 Input Parameters:
1951 + mat - a dense matrix
1952 - array - pointer to the data
1953
1954 Level: intermediate
1955
1956 .seealso: MatDenseGetArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite()
1957 @*/
MatDenseRestoreArrayRead(Mat A,const PetscScalar ** array)1958 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array)
1959 {
1960 PetscErrorCode ierr;
1961
1962 PetscFunctionBegin;
1963 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
1964 PetscValidPointer(array,2);
1965 ierr = PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr);
1966 PetscFunctionReturn(0);
1967 }
1968
1969 /*@C
1970 MatDenseGetArrayWrite - gives write-only access to the array where the data for a dense matrix is stored
1971
1972 Not Collective
1973
1974 Input Parameter:
1975 . mat - a dense matrix
1976
1977 Output Parameter:
1978 . array - pointer to the data
1979
1980 Level: intermediate
1981
1982 .seealso: MatDenseRestoreArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
1983 @*/
MatDenseGetArrayWrite(Mat A,PetscScalar ** array)1984 PetscErrorCode MatDenseGetArrayWrite(Mat A,PetscScalar **array)
1985 {
1986 PetscErrorCode ierr;
1987
1988 PetscFunctionBegin;
1989 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
1990 PetscValidPointer(array,2);
1991 ierr = PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr);
1992 PetscFunctionReturn(0);
1993 }
1994
1995 /*@C
1996 MatDenseRestoreArrayWrite - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayWrite()
1997
1998 Not Collective
1999
2000 Input Parameters:
2001 + mat - a dense matrix
2002 - array - pointer to the data
2003
2004 Level: intermediate
2005
2006 .seealso: MatDenseGetArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead()
2007 @*/
MatDenseRestoreArrayWrite(Mat A,PetscScalar ** array)2008 PetscErrorCode MatDenseRestoreArrayWrite(Mat A,PetscScalar **array)
2009 {
2010 PetscErrorCode ierr;
2011
2012 PetscFunctionBegin;
2013 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
2014 PetscValidPointer(array,2);
2015 ierr = PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr);
2016 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr);
2017 #if defined(PETSC_HAVE_CUDA)
2018 A->offloadmask = PETSC_OFFLOAD_CPU;
2019 #endif
2020 PetscFunctionReturn(0);
2021 }
2022
MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat * B)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;
2026 PetscErrorCode ierr;
2027 PetscInt i,j,nrows,ncols,blda;
2028 const PetscInt *irow,*icol;
2029 PetscScalar *av,*bv,*v = mat->v;
2030 Mat newmat;
2031
2032 PetscFunctionBegin;
2033 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr);
2034 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr);
2035 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr);
2036 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr);
2037
2038 /* Check submatrixcall */
2039 if (scall == MAT_REUSE_MATRIX) {
2040 PetscInt n_cols,n_rows;
2041 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr);
2042 if (n_rows != nrows || n_cols != ncols) {
2043 /* resize the result matrix to match number of requested rows/columns */
2044 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr);
2045 }
2046 newmat = *B;
2047 } else {
2048 /* Create and fill new matrix */
2049 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr);
2050 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr);
2051 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr);
2052 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr);
2053 }
2054
2055 /* Now extract the data pointers and do the copy,column at a time */
2056 ierr = MatDenseGetArray(newmat,&bv);CHKERRQ(ierr);
2057 ierr = MatDenseGetLDA(newmat,&blda);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArray(newmat,&bv);CHKERRQ(ierr);
2064
2065 /* Assemble the matrices so that the correct flags are set */
2066 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2067 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2068
2069 /* Free work space */
2070 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr);
2071 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr);
2072 *B = newmat;
2073 PetscFunctionReturn(0);
2074 }
2075
MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat * B[])2076 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[])
2077 {
2078 PetscErrorCode ierr;
2079 PetscInt i;
2080
2081 PetscFunctionBegin;
2082 if (scall == MAT_INITIAL_MATRIX) {
2083 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr);
2084 }
2085
2086 for (i=0; i<n; i++) {
2087 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr);
2088 }
2089 PetscFunctionReturn(0);
2090 }
2091
MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)2092 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode)
2093 {
2094 PetscFunctionBegin;
2095 PetscFunctionReturn(0);
2096 }
2097
MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)2098 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode)
2099 {
2100 PetscFunctionBegin;
2101 PetscFunctionReturn(0);
2102 }
2103
MatCopy_SeqDense(Mat A,Mat B,MatStructure str)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;
2111
2112 PetscFunctionBegin;
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 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr);
2116 PetscFunctionReturn(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 ierr = MatDenseGetArrayRead(A,&va);CHKERRQ(ierr);
2120 ierr = MatDenseGetArray(B,&vb);CHKERRQ(ierr);
2121 if (lda1>m || lda2>m) {
2122 for (j=0; j<n; j++) {
2123 ierr = PetscArraycpy(vb+j*lda2,va+j*lda1,m);CHKERRQ(ierr);
2124 }
2125 } else {
2126 ierr = PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);CHKERRQ(ierr);
2127 }
2128 ierr = MatDenseRestoreArray(B,&vb);CHKERRQ(ierr);
2129 ierr = MatDenseRestoreArrayRead(A,&va);CHKERRQ(ierr);
2130 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2131 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2132 PetscFunctionReturn(0);
2133 }
2134
MatSetUp_SeqDense(Mat A)2135 static PetscErrorCode MatSetUp_SeqDense(Mat A)
2136 {
2137 PetscErrorCode ierr;
2138
2139 PetscFunctionBegin;
2140 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr);
2141 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr);
2142 if (!A->preallocated) {
2143 ierr = MatSeqDenseSetPreallocation(A,NULL);CHKERRQ(ierr);
2144 }
2145 PetscFunctionReturn(0);
2146 }
2147
MatConjugate_SeqDense(Mat A)2148 static PetscErrorCode MatConjugate_SeqDense(Mat A)
2149 {
2150 PetscInt i,nz = A->rmap->n*A->cmap->n;
2151 PetscScalar *aa;
2152 PetscErrorCode ierr;
2153
2154 PetscFunctionBegin;
2155 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr);
2156 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]);
2157 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr);
2158 PetscFunctionReturn(0);
2159 }
2160
MatRealPart_SeqDense(Mat A)2161 static PetscErrorCode MatRealPart_SeqDense(Mat A)
2162 {
2163 PetscInt i,nz = A->rmap->n*A->cmap->n;
2164 PetscScalar *aa;
2165 PetscErrorCode ierr;
2166
2167 PetscFunctionBegin;
2168 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr);
2169 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
2170 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr);
2171 PetscFunctionReturn(0);
2172 }
2173
MatImaginaryPart_SeqDense(Mat A)2174 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2175 {
2176 PetscInt i,nz = A->rmap->n*A->cmap->n;
2177 PetscScalar *aa;
2178 PetscErrorCode ierr;
2179
2180 PetscFunctionBegin;
2181 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr);
2182 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2183 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr);
2184 PetscFunctionReturn(0);
2185 }
2186
2187 /* ----------------------------------------------------------------*/
MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)2188 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2189 {
2190 PetscErrorCode ierr;
2191 PetscInt m=A->rmap->n,n=B->cmap->n;
2192 PetscBool cisdense;
2193
2194 PetscFunctionBegin;
2195 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr);
2196 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr);
2197 if (!cisdense) {
2198 PetscBool flg;
2199
2200 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr);
2201 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr);
2202 }
2203 ierr = MatSetUp(C);CHKERRQ(ierr);
2204 PetscFunctionReturn(0);
2205 }
2206
MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)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;
2215
2216 PetscFunctionBegin;
2217 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr);
2218 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr);
2219 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr);
2220 if (!m || !n || !k) PetscFunctionReturn(0);
2221 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr);
2222 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr);
2223 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr);
2224 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2225 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr);
2226 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr);
2227 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr);
2228 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr);
2229 PetscFunctionReturn(0);
2230 }
2231
MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)2232 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2233 {
2234 PetscErrorCode ierr;
2235 PetscInt m=A->rmap->n,n=B->rmap->n;
2236 PetscBool cisdense;
2237
2238 PetscFunctionBegin;
2239 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr);
2240 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr);
2241 if (!cisdense) {
2242 PetscBool flg;
2243
2244 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr);
2245 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr);
2246 }
2247 ierr = MatSetUp(C);CHKERRQ(ierr);
2248 PetscFunctionReturn(0);
2249 }
2250
MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)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;
2261
2262 PetscFunctionBegin;
2263 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr);
2264 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr);
2265 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr);
2266 if (!m || !n || !k) PetscFunctionReturn(0);
2267 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr);
2268 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr);
2269 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr);
2270 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2271 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr);
2272 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr);
2273 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr);
2274 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr);
2275 PetscFunctionReturn(0);
2276 }
2277
MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)2278 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C)
2279 {
2280 PetscErrorCode ierr;
2281 PetscInt m=A->cmap->n,n=B->cmap->n;
2282 PetscBool cisdense;
2283
2284 PetscFunctionBegin;
2285 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr);
2286 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr);
2287 if (!cisdense) {
2288 PetscBool flg;
2289
2290 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr);
2291 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr);
2292 }
2293 ierr = MatSetUp(C);CHKERRQ(ierr);
2294 PetscFunctionReturn(0);
2295 }
2296
MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C)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;
2307
2308 PetscFunctionBegin;
2309 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr);
2310 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr);
2311 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr);
2312 if (!m || !n || !k) PetscFunctionReturn(0);
2313 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr);
2314 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr);
2315 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr);
2316 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda));
2317 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr);
2318 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr);
2319 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr);
2320 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr);
2321 PetscFunctionReturn(0);
2322 }
2323
2324 /* ----------------------------------------------- */
MatProductSetFromOptions_SeqDense_AB(Mat C)2325 static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2326 {
2327 PetscFunctionBegin;
2328 C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2329 C->ops->productsymbolic = MatProductSymbolic_AB;
2330 PetscFunctionReturn(0);
2331 }
2332
MatProductSetFromOptions_SeqDense_AtB(Mat C)2333 static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2334 {
2335 PetscFunctionBegin;
2336 C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2337 C->ops->productsymbolic = MatProductSymbolic_AtB;
2338 PetscFunctionReturn(0);
2339 }
2340
MatProductSetFromOptions_SeqDense_ABt(Mat C)2341 static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2342 {
2343 PetscFunctionBegin;
2344 C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2345 C->ops->productsymbolic = MatProductSymbolic_ABt;
2346 PetscFunctionReturn(0);
2347 }
2348
MatProductSetFromOptions_SeqDense(Mat C)2349 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2350 {
2351 PetscErrorCode ierr;
2352 Mat_Product *product = C->product;
2353
2354 PetscFunctionBegin;
2355 switch (product->type) {
2356 case MATPRODUCT_AB:
2357 ierr = MatProductSetFromOptions_SeqDense_AB(C);CHKERRQ(ierr);
2358 break;
2359 case MATPRODUCT_AtB:
2360 ierr = MatProductSetFromOptions_SeqDense_AtB(C);CHKERRQ(ierr);
2361 break;
2362 case MATPRODUCT_ABt:
2363 ierr = MatProductSetFromOptions_SeqDense_ABt(C);CHKERRQ(ierr);
2364 break;
2365 default:
2366 break;
2367 }
2368 PetscFunctionReturn(0);
2369 }
2370 /* ----------------------------------------------- */
2371
MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[])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;
2379
2380 PetscFunctionBegin;
2381 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2382 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
2383 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr);
2384 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr);
2393 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
2394 PetscFunctionReturn(0);
2395 }
2396
MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[])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;
2405
2406 PetscFunctionBegin;
2407 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2408 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
2409 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr);
2410 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr);
2420 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
2421 PetscFunctionReturn(0);
2422 }
2423
MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[])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;
2431
2432 PetscFunctionBegin;
2433 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2434 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr);
2435 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
2436 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr);
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 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
2445 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr);
2446 PetscFunctionReturn(0);
2447 }
2448
MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col)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;
2455
2456 PetscFunctionBegin;
2457 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2458 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr);
2459 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
2460 ierr = PetscArraycpy(x,aa+col*a->lda,A->rmap->n);CHKERRQ(ierr);
2461 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
2462 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr);
2463 PetscFunctionReturn(0);
2464 }
2465
MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal * norms)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;
2471
2472 PetscFunctionBegin;
2473 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);
2474 ierr = PetscArrayzero(norms,n);CHKERRQ(ierr);
2475 ierr = MatDenseGetArrayRead(A,&a);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArrayRead(A,&a);CHKERRQ(ierr);
2499 if (type == NORM_2) {
2500 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]);
2501 }
2502 PetscFunctionReturn(0);
2503 }
2504
MatSetRandom_SeqDense(Mat x,PetscRandom rctx)2505 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx)
2506 {
2507 PetscErrorCode ierr;
2508 PetscScalar *a;
2509 PetscInt lda,m,n,i,j;
2510
2511 PetscFunctionBegin;
2512 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr);
2513 ierr = MatDenseGetLDA(x,&lda);CHKERRQ(ierr);
2514 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr);
2515 for (j=0; j<n; j++) {
2516 for (i=0; i<m; i++) {
2517 ierr = PetscRandomGetValue(rctx,a+j*lda+i);CHKERRQ(ierr);
2518 }
2519 }
2520 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr);
2521 PetscFunctionReturn(0);
2522 }
2523
MatMissingDiagonal_SeqDense(Mat A,PetscBool * missing,PetscInt * d)2524 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d)
2525 {
2526 PetscFunctionBegin;
2527 *missing = PETSC_FALSE;
2528 PetscFunctionReturn(0);
2529 }
2530
2531 /* vals is not const */
MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar ** vals)2532 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals)
2533 {
2534 PetscErrorCode ierr;
2535 Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2536 PetscScalar *v;
2537
2538 PetscFunctionBegin;
2539 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2540 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr);
2541 *vals = v+col*a->lda;
2542 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr);
2543 PetscFunctionReturn(0);
2544 }
2545
MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar ** vals)2546 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals)
2547 {
2548 PetscFunctionBegin;
2549 *vals = NULL; /* user cannot accidently use the array later */
2550 PetscFunctionReturn(0);
2551 }
2552
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 };
2703
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.
2708
2709 Collective
2710
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.
2717
2718 Output Parameter:
2719 . A - the matrix
2720
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.
2725
2726 Level: intermediate
2727
2728 .seealso: MatCreate(), MatCreateDense(), MatSetValues()
2729 @*/
MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar * data,Mat * A)2730 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A)
2731 {
2732 PetscErrorCode ierr;
2733
2734 PetscFunctionBegin;
2735 ierr = MatCreate(comm,A);CHKERRQ(ierr);
2736 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr);
2737 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr);
2738 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr);
2739 PetscFunctionReturn(0);
2740 }
2741
2742 /*@C
2743 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements
2744
2745 Collective
2746
2747 Input Parameters:
2748 + B - the matrix
2749 - data - the array (or NULL)
2750
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.
2755
2756 Level: intermediate
2757
2758 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatDenseSetLDA()
2759
2760 @*/
MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])2761 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[])
2762 {
2763 PetscErrorCode ierr;
2764
2765 PetscFunctionBegin;
2766 PetscValidHeaderSpecific(B,MAT_CLASSID,1);
2767 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr);
2768 PetscFunctionReturn(0);
2769 }
2770
MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar * data)2771 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data)
2772 {
2773 Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2774 PetscErrorCode ierr;
2775
2776 PetscFunctionBegin;
2777 if (b->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first");
2778 B->preallocated = PETSC_TRUE;
2779
2780 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr);
2781 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr);
2782
2783 if (b->lda <= 0) b->lda = B->rmap->n;
2784
2785 ierr = PetscIntMultError(b->lda,B->cmap->n,NULL);CHKERRQ(ierr);
2786 if (!data) { /* petsc-allocated storage */
2787 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); }
2788 ierr = PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v);CHKERRQ(ierr);
2789 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr);
2790
2791 b->user_alloc = PETSC_FALSE;
2792 } else { /* user-allocated storage */
2793 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); }
2794 b->v = data;
2795 b->user_alloc = PETSC_TRUE;
2796 }
2797 B->assembled = PETSC_TRUE;
2798 PetscFunctionReturn(0);
2799 }
2800
2801 #if defined(PETSC_HAVE_ELEMENTAL)
MatConvert_SeqDense_Elemental(Mat A,MatType newtype,MatReuse reuse,Mat * newmat)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;
2809
2810 PetscFunctionBegin;
2811 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr);
2812 ierr = MatDenseGetArrayRead(A,&array);CHKERRQ(ierr);
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 ierr = MatDenseRestoreArrayRead(A,&array);CHKERRQ(ierr);
2825
2826 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr);
2827 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr);
2828 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr);
2829 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr);
2830
2831 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
2832 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr);
2833 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2834 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2835 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr);
2836
2837 if (reuse == MAT_INPLACE_MATRIX) {
2838 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr);
2839 } else {
2840 *newmat = mat_elemental;
2841 }
2842 PetscFunctionReturn(0);
2843 }
2844 #endif
2845
MatDenseSetLDA_SeqDense(Mat B,PetscInt lda)2846 PetscErrorCode MatDenseSetLDA_SeqDense(Mat B,PetscInt lda)
2847 {
2848 Mat_SeqDense *b = (Mat_SeqDense*)B->data;
2849 PetscBool data;
2850
2851 PetscFunctionBegin;
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 PetscFunctionReturn(0);
2857 }
2858
MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat * outmat)2859 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
2860 {
2861 PetscErrorCode ierr;
2862 PetscMPIInt size;
2863
2864 PetscFunctionBegin;
2865 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
2866 if (size == 1) {
2867 if (scall == MAT_INITIAL_MATRIX) {
2868 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr);
2869 } else {
2870 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
2871 }
2872 } else {
2873 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr);
2874 }
2875 PetscFunctionReturn(0);
2876 }
2877
MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec * v)2878 PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
2879 {
2880 Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2881 PetscErrorCode ierr;
2882
2883 PetscFunctionBegin;
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 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr);
2888 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr);
2889 }
2890 a->vecinuse = col + 1;
2891 ierr = MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr);
2892 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr);
2893 *v = a->cvec;
2894 PetscFunctionReturn(0);
2895 }
2896
MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec * v)2897 PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v)
2898 {
2899 Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2900 PetscErrorCode ierr;
2901
2902 PetscFunctionBegin;
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 ierr = MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr);
2907 ierr = VecResetArray(a->cvec);CHKERRQ(ierr);
2908 *v = NULL;
2909 PetscFunctionReturn(0);
2910 }
2911
MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec * v)2912 PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
2913 {
2914 Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2915 PetscErrorCode ierr;
2916
2917 PetscFunctionBegin;
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 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr);
2922 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr);
2923 }
2924 a->vecinuse = col + 1;
2925 ierr = MatDenseGetArrayRead(A,&a->ptrinuse);CHKERRQ(ierr);
2926 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr);
2927 ierr = VecLockReadPush(a->cvec);CHKERRQ(ierr);
2928 *v = a->cvec;
2929 PetscFunctionReturn(0);
2930 }
2931
MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec * v)2932 PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v)
2933 {
2934 Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2935 PetscErrorCode ierr;
2936
2937 PetscFunctionBegin;
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 ierr = MatDenseRestoreArrayRead(A,&a->ptrinuse);CHKERRQ(ierr);
2942 ierr = VecLockReadPop(a->cvec);CHKERRQ(ierr);
2943 ierr = VecResetArray(a->cvec);CHKERRQ(ierr);
2944 *v = NULL;
2945 PetscFunctionReturn(0);
2946 }
2947
MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec * v)2948 PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
2949 {
2950 Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2951 PetscErrorCode ierr;
2952
2953 PetscFunctionBegin;
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 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr);
2958 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr);
2959 }
2960 a->vecinuse = col + 1;
2961 ierr = MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr);
2962 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr);
2963 *v = a->cvec;
2964 PetscFunctionReturn(0);
2965 }
2966
MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec * v)2967 PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v)
2968 {
2969 Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2970 PetscErrorCode ierr;
2971
2972 PetscFunctionBegin;
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 ierr = MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr);
2977 ierr = VecResetArray(a->cvec);CHKERRQ(ierr);
2978 *v = NULL;
2979 PetscFunctionReturn(0);
2980 }
2981
MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat * v)2982 PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
2983 {
2984 Mat_SeqDense *a = (Mat_SeqDense*)A->data;
2985 PetscErrorCode ierr;
2986
2987 PetscFunctionBegin;
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 ierr = MatDestroy(&a->cmat);CHKERRQ(ierr);
2992 }
2993 if (!a->cmat) {
2994 ierr = 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);CHKERRQ(ierr);
2995 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat);CHKERRQ(ierr);
2996 } else {
2997 ierr = MatDensePlaceArray(a->cmat,a->v + (size_t)cbegin * (size_t)a->lda);CHKERRQ(ierr);
2998 }
2999 ierr = MatDenseSetLDA(a->cmat,a->lda);CHKERRQ(ierr);
3000 a->matinuse = cbegin + 1;
3001 *v = a->cmat;
3002 PetscFunctionReturn(0);
3003 }
3004
MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat * v)3005 PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v)
3006 {
3007 Mat_SeqDense *a = (Mat_SeqDense*)A->data;
3008 PetscErrorCode ierr;
3009
3010 PetscFunctionBegin;
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 ierr = MatDenseResetArray(a->cmat);CHKERRQ(ierr);
3016 *v = NULL;
3017 PetscFunctionReturn(0);
3018 }
3019
3020 /*MC
3021 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
3022
3023 Options Database Keys:
3024 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions()
3025
3026 Level: beginner
3027
3028 .seealso: MatCreateSeqDense()
3029
3030 M*/
MatCreate_SeqDense(Mat B)3031 PetscErrorCode MatCreate_SeqDense(Mat B)
3032 {
3033 Mat_SeqDense *b;
3034 PetscErrorCode ierr;
3035 PetscMPIInt size;
3036
3037 PetscFunctionBegin;
3038 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr);
3039 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");
3040
3041 ierr = PetscNewLog(B,&b);CHKERRQ(ierr);
3042 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr);
3043 B->data = (void*)b;
3044
3045 b->roworiented = PETSC_TRUE;
3046
3047 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);CHKERRQ(ierr);
3048 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense);CHKERRQ(ierr);
3049 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr);
3050 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr);
3051 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr);
3052 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr);
3053 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense);CHKERRQ(ierr);
3054 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr);
3055 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr);
3056 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr);
3057 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr);
3058 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr);
3059 #if defined(PETSC_HAVE_ELEMENTAL)
3060 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr);
3061 #endif
3062 #if defined(PETSC_HAVE_SCALAPACK)
3063 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_scalapack_C",MatConvert_Dense_ScaLAPACK);CHKERRQ(ierr);
3064 #endif
3065 #if defined(PETSC_HAVE_CUDA)
3066 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);CHKERRQ(ierr);
3067 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr);
3068 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr);
3069 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr);
3070 #endif
3071 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr);
3072 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);CHKERRQ(ierr);
3073 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr);
3074 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr);
3075 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr);
3076
3077 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr);
3078 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr);
3079 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);CHKERRQ(ierr);
3080 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);CHKERRQ(ierr);
3081 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);CHKERRQ(ierr);
3082 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);CHKERRQ(ierr);
3083 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);CHKERRQ(ierr);
3084 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);CHKERRQ(ierr);
3085 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense);CHKERRQ(ierr);
3086 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense);CHKERRQ(ierr);
3087 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr);
3088 PetscFunctionReturn(0);
3089 }
3090
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.
3093
3094 Not Collective
3095
3096 Input Parameters:
3097 + mat - a MATSEQDENSE or MATMPIDENSE matrix
3098 - col - column index
3099
3100 Output Parameter:
3101 . vals - pointer to the data
3102
3103 Level: intermediate
3104
3105 .seealso: MatDenseRestoreColumn()
3106 @*/
MatDenseGetColumn(Mat A,PetscInt col,PetscScalar ** vals)3107 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals)
3108 {
3109 PetscErrorCode ierr;
3110
3111 PetscFunctionBegin;
3112 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
3113 PetscValidLogicalCollectiveInt(A,col,2);
3114 PetscValidPointer(vals,3);
3115 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr);
3116 PetscFunctionReturn(0);
3117 }
3118
3119 /*@C
3120 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn().
3121
3122 Not Collective
3123
3124 Input Parameter:
3125 . mat - a MATSEQDENSE or MATMPIDENSE matrix
3126
3127 Output Parameter:
3128 . vals - pointer to the data
3129
3130 Level: intermediate
3131
3132 .seealso: MatDenseGetColumn()
3133 @*/
MatDenseRestoreColumn(Mat A,PetscScalar ** vals)3134 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals)
3135 {
3136 PetscErrorCode ierr;
3137
3138 PetscFunctionBegin;
3139 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
3140 PetscValidPointer(vals,2);
3141 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr);
3142 PetscFunctionReturn(0);
3143 }
3144
3145 /*@C
3146 MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec.
3147
3148 Collective
3149
3150 Input Parameters:
3151 + mat - the Mat object
3152 - col - the column index
3153
3154 Output Parameter:
3155 . v - the vector
3156
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.
3160
3161 Level: intermediate
3162
3163 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3164 @*/
MatDenseGetColumnVec(Mat A,PetscInt col,Vec * v)3165 PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v)
3166 {
3167 PetscErrorCode ierr;
3168
3169 PetscFunctionBegin;
3170 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
3171 PetscValidType(A,1);
3172 PetscValidLogicalCollectiveInt(A,col,2);
3173 PetscValidPointer(v,3);
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 ierr = PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr);
3177 PetscFunctionReturn(0);
3178 }
3179
3180 /*@C
3181 MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec().
3182
3183 Collective
3184
3185 Input Parameters:
3186 + mat - the Mat object
3187 . col - the column index
3188 - v - the Vec object
3189
3190 Level: intermediate
3191
3192 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3193 @*/
MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec * v)3194 PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v)
3195 {
3196 PetscErrorCode ierr;
3197
3198 PetscFunctionBegin;
3199 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
3200 PetscValidType(A,1);
3201 PetscValidLogicalCollectiveInt(A,col,2);
3202 PetscValidPointer(v,3);
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 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr);
3206 PetscFunctionReturn(0);
3207 }
3208
3209 /*@C
3210 MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec.
3211
3212 Collective
3213
3214 Input Parameters:
3215 + mat - the Mat object
3216 - col - the column index
3217
3218 Output Parameter:
3219 . v - the vector
3220
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.
3225
3226 Level: intermediate
3227
3228 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3229 @*/
MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec * v)3230 PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v)
3231 {
3232 PetscErrorCode ierr;
3233
3234 PetscFunctionBegin;
3235 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
3236 PetscValidType(A,1);
3237 PetscValidLogicalCollectiveInt(A,col,2);
3238 PetscValidPointer(v,3);
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 ierr = PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr);
3242 PetscFunctionReturn(0);
3243 }
3244
3245 /*@C
3246 MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead().
3247
3248 Collective
3249
3250 Input Parameters:
3251 + mat - the Mat object
3252 . col - the column index
3253 - v - the Vec object
3254
3255 Level: intermediate
3256
3257 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite()
3258 @*/
MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec * v)3259 PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v)
3260 {
3261 PetscErrorCode ierr;
3262
3263 PetscFunctionBegin;
3264 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
3265 PetscValidType(A,1);
3266 PetscValidLogicalCollectiveInt(A,col,2);
3267 PetscValidPointer(v,3);
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 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr);
3271 PetscFunctionReturn(0);
3272 }
3273
3274 /*@C
3275 MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec.
3276
3277 Collective
3278
3279 Input Parameters:
3280 + mat - the Mat object
3281 - col - the column index
3282
3283 Output Parameter:
3284 . v - the vector
3285
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.
3289
3290 Level: intermediate
3291
3292 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite()
3293 @*/
MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec * v)3294 PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v)
3295 {
3296 PetscErrorCode ierr;
3297
3298 PetscFunctionBegin;
3299 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
3300 PetscValidType(A,1);
3301 PetscValidLogicalCollectiveInt(A,col,2);
3302 PetscValidPointer(v,3);
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 ierr = PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr);
3306 PetscFunctionReturn(0);
3307 }
3308
3309 /*@C
3310 MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite().
3311
3312 Collective
3313
3314 Input Parameters:
3315 + mat - the Mat object
3316 . col - the column index
3317 - v - the Vec object
3318
3319 Level: intermediate
3320
3321 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead()
3322 @*/
MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec * v)3323 PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v)
3324 {
3325 PetscErrorCode ierr;
3326
3327 PetscFunctionBegin;
3328 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
3329 PetscValidType(A,1);
3330 PetscValidLogicalCollectiveInt(A,col,2);
3331 PetscValidPointer(v,3);
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 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr);
3335 PetscFunctionReturn(0);
3336 }
3337
3338 /*@C
3339 MatDenseGetSubMatrix - Gives access to a block of columns of a dense matrix, represented as a Mat.
3340
3341 Collective
3342
3343 Input Parameters:
3344 + mat - the Mat object
3345 . cbegin - the first index in the block
3346 - cend - the last index in the block
3347
3348 Output Parameter:
3349 . v - the matrix
3350
3351 Notes:
3352 The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed.
3353
3354 Level: intermediate
3355
3356 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseRestoreSubMatrix()
3357 @*/
MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat * v)3358 PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat *v)
3359 {
3360 PetscErrorCode ierr;
3361
3362 PetscFunctionBegin;
3363 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
3364 PetscValidType(A,1);
3365 PetscValidLogicalCollectiveInt(A,cbegin,2);
3366 PetscValidLogicalCollectiveInt(A,cend,3);
3367 PetscValidPointer(v,4);
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 ierr = PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));CHKERRQ(ierr);
3372 PetscFunctionReturn(0);
3373 }
3374
3375 /*@C
3376 MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix().
3377
3378 Collective
3379
3380 Input Parameters:
3381 + mat - the Mat object
3382 - v - the Mat object
3383
3384 Level: intermediate
3385
3386 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseGetSubMatrix()
3387 @*/
MatDenseRestoreSubMatrix(Mat A,Mat * v)3388 PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v)
3389 {
3390 PetscErrorCode ierr;
3391
3392 PetscFunctionBegin;
3393 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
3394 PetscValidType(A,1);
3395 PetscValidPointer(v,2);
3396 ierr = PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v));CHKERRQ(ierr);
3397 PetscFunctionReturn(0);
3398 }
3399