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