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