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