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