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 case MATPRODUCT_RARt: 2309 ierr = MatProductSetFromOptions_SeqDense_PtAP(C);CHKERRQ(ierr); 2310 break; 2311 default: SETERRQ1(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatProduct type %s is not supported for SeqDense and SeqDense matrices",MatProductTypes[product->type]); 2312 } 2313 PetscFunctionReturn(0); 2314 } 2315 /* ----------------------------------------------- */ 2316 2317 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2318 { 2319 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2320 PetscErrorCode ierr; 2321 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2322 PetscScalar *x; 2323 const PetscScalar *aa; 2324 2325 PetscFunctionBegin; 2326 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2327 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2328 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2329 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2330 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2331 for (i=0; i<m; i++) { 2332 x[i] = aa[i]; if (idx) idx[i] = 0; 2333 for (j=1; j<n; j++) { 2334 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2335 } 2336 } 2337 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2338 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2339 PetscFunctionReturn(0); 2340 } 2341 2342 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2343 { 2344 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2345 PetscErrorCode ierr; 2346 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2347 PetscScalar *x; 2348 PetscReal atmp; 2349 const PetscScalar *aa; 2350 2351 PetscFunctionBegin; 2352 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2353 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2354 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2355 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2356 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2357 for (i=0; i<m; i++) { 2358 x[i] = PetscAbsScalar(aa[i]); 2359 for (j=1; j<n; j++) { 2360 atmp = PetscAbsScalar(aa[i+a->lda*j]); 2361 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2362 } 2363 } 2364 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2365 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2366 PetscFunctionReturn(0); 2367 } 2368 2369 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2370 { 2371 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2372 PetscErrorCode ierr; 2373 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2374 PetscScalar *x; 2375 const PetscScalar *aa; 2376 2377 PetscFunctionBegin; 2378 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2379 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2380 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2381 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2382 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2383 for (i=0; i<m; i++) { 2384 x[i] = aa[i]; if (idx) idx[i] = 0; 2385 for (j=1; j<n; j++) { 2386 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2387 } 2388 } 2389 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2390 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2391 PetscFunctionReturn(0); 2392 } 2393 2394 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2395 { 2396 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2397 PetscErrorCode ierr; 2398 PetscScalar *x; 2399 const PetscScalar *aa; 2400 2401 PetscFunctionBegin; 2402 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2403 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2404 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2405 ierr = PetscArraycpy(x,aa+col*a->lda,A->rmap->n);CHKERRQ(ierr); 2406 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2407 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2408 PetscFunctionReturn(0); 2409 } 2410 2411 PETSC_INTERN PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2412 { 2413 PetscErrorCode ierr; 2414 PetscInt i,j,m,n; 2415 const PetscScalar *a; 2416 2417 PetscFunctionBegin; 2418 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2419 ierr = PetscArrayzero(norms,n);CHKERRQ(ierr); 2420 ierr = MatDenseGetArrayRead(A,&a);CHKERRQ(ierr); 2421 if (type == NORM_2) { 2422 for (i=0; i<n; i++) { 2423 for (j=0; j<m; j++) { 2424 norms[i] += PetscAbsScalar(a[j]*a[j]); 2425 } 2426 a += m; 2427 } 2428 } else if (type == NORM_1) { 2429 for (i=0; i<n; i++) { 2430 for (j=0; j<m; j++) { 2431 norms[i] += PetscAbsScalar(a[j]); 2432 } 2433 a += m; 2434 } 2435 } else if (type == NORM_INFINITY) { 2436 for (i=0; i<n; i++) { 2437 for (j=0; j<m; j++) { 2438 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2439 } 2440 a += m; 2441 } 2442 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2443 ierr = MatDenseRestoreArrayRead(A,&a);CHKERRQ(ierr); 2444 if (type == NORM_2) { 2445 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2446 } 2447 PetscFunctionReturn(0); 2448 } 2449 2450 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2451 { 2452 PetscErrorCode ierr; 2453 PetscScalar *a; 2454 PetscInt lda,m,n,i,j; 2455 2456 PetscFunctionBegin; 2457 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2458 ierr = MatDenseGetLDA(x,&lda);CHKERRQ(ierr); 2459 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2460 for (j=0; j<n; j++) { 2461 for (i=0; i<m; i++) { 2462 ierr = PetscRandomGetValue(rctx,a+j*lda+i);CHKERRQ(ierr); 2463 } 2464 } 2465 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2466 PetscFunctionReturn(0); 2467 } 2468 2469 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2470 { 2471 PetscFunctionBegin; 2472 *missing = PETSC_FALSE; 2473 PetscFunctionReturn(0); 2474 } 2475 2476 /* vals is not const */ 2477 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2478 { 2479 PetscErrorCode ierr; 2480 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2481 PetscScalar *v; 2482 2483 PetscFunctionBegin; 2484 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2485 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 2486 *vals = v+col*a->lda; 2487 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 2488 PetscFunctionReturn(0); 2489 } 2490 2491 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2492 { 2493 PetscFunctionBegin; 2494 *vals = 0; /* user cannot accidently use the array later */ 2495 PetscFunctionReturn(0); 2496 } 2497 2498 /* -------------------------------------------------------------------*/ 2499 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2500 MatGetRow_SeqDense, 2501 MatRestoreRow_SeqDense, 2502 MatMult_SeqDense, 2503 /* 4*/ MatMultAdd_SeqDense, 2504 MatMultTranspose_SeqDense, 2505 MatMultTransposeAdd_SeqDense, 2506 0, 2507 0, 2508 0, 2509 /* 10*/ 0, 2510 MatLUFactor_SeqDense, 2511 MatCholeskyFactor_SeqDense, 2512 MatSOR_SeqDense, 2513 MatTranspose_SeqDense, 2514 /* 15*/ MatGetInfo_SeqDense, 2515 MatEqual_SeqDense, 2516 MatGetDiagonal_SeqDense, 2517 MatDiagonalScale_SeqDense, 2518 MatNorm_SeqDense, 2519 /* 20*/ MatAssemblyBegin_SeqDense, 2520 MatAssemblyEnd_SeqDense, 2521 MatSetOption_SeqDense, 2522 MatZeroEntries_SeqDense, 2523 /* 24*/ MatZeroRows_SeqDense, 2524 0, 2525 0, 2526 0, 2527 0, 2528 /* 29*/ MatSetUp_SeqDense, 2529 0, 2530 0, 2531 0, 2532 0, 2533 /* 34*/ MatDuplicate_SeqDense, 2534 0, 2535 0, 2536 0, 2537 0, 2538 /* 39*/ MatAXPY_SeqDense, 2539 MatCreateSubMatrices_SeqDense, 2540 0, 2541 MatGetValues_SeqDense, 2542 MatCopy_SeqDense, 2543 /* 44*/ MatGetRowMax_SeqDense, 2544 MatScale_SeqDense, 2545 MatShift_Basic, 2546 0, 2547 MatZeroRowsColumns_SeqDense, 2548 /* 49*/ MatSetRandom_SeqDense, 2549 0, 2550 0, 2551 0, 2552 0, 2553 /* 54*/ 0, 2554 0, 2555 0, 2556 0, 2557 0, 2558 /* 59*/ 0, 2559 MatDestroy_SeqDense, 2560 MatView_SeqDense, 2561 0, 2562 0, 2563 /* 64*/ 0, 2564 0, 2565 0, 2566 0, 2567 0, 2568 /* 69*/ MatGetRowMaxAbs_SeqDense, 2569 0, 2570 0, 2571 0, 2572 0, 2573 /* 74*/ 0, 2574 0, 2575 0, 2576 0, 2577 0, 2578 /* 79*/ 0, 2579 0, 2580 0, 2581 0, 2582 /* 83*/ MatLoad_SeqDense, 2583 MatIsSymmetric_SeqDense, 2584 MatIsHermitian_SeqDense, 2585 0, 2586 0, 2587 0, 2588 /* 89*/ 0, 2589 0, 2590 MatMatMultNumeric_SeqDense_SeqDense, 2591 0, 2592 0, 2593 /* 94*/ 0, 2594 0, 2595 0, 2596 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2597 0, 2598 /* 99*/ MatProductSetFromOptions_SeqDense, 2599 0, 2600 0, 2601 MatConjugate_SeqDense, 2602 0, 2603 /*104*/ 0, 2604 MatRealPart_SeqDense, 2605 MatImaginaryPart_SeqDense, 2606 0, 2607 0, 2608 /*109*/ 0, 2609 0, 2610 MatGetRowMin_SeqDense, 2611 MatGetColumnVector_SeqDense, 2612 MatMissingDiagonal_SeqDense, 2613 /*114*/ 0, 2614 0, 2615 0, 2616 0, 2617 0, 2618 /*119*/ 0, 2619 0, 2620 0, 2621 0, 2622 0, 2623 /*124*/ 0, 2624 MatGetColumnNorms_SeqDense, 2625 0, 2626 0, 2627 0, 2628 /*129*/ 0, 2629 0, 2630 0, 2631 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2632 0, 2633 /*134*/ 0, 2634 0, 2635 0, 2636 0, 2637 0, 2638 /*139*/ 0, 2639 0, 2640 0, 2641 0, 2642 0, 2643 MatCreateMPIMatConcatenateSeqMat_SeqDense, 2644 /*145*/ 0, 2645 0, 2646 0 2647 }; 2648 2649 /*@C 2650 MatCreateSeqDense - Creates a sequential dense matrix that 2651 is stored in column major order (the usual Fortran 77 manner). Many 2652 of the matrix operations use the BLAS and LAPACK routines. 2653 2654 Collective 2655 2656 Input Parameters: 2657 + comm - MPI communicator, set to PETSC_COMM_SELF 2658 . m - number of rows 2659 . n - number of columns 2660 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2661 to control all matrix memory allocation. 2662 2663 Output Parameter: 2664 . A - the matrix 2665 2666 Notes: 2667 The data input variable is intended primarily for Fortran programmers 2668 who wish to allocate their own matrix memory space. Most users should 2669 set data=NULL. 2670 2671 Level: intermediate 2672 2673 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2674 @*/ 2675 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2676 { 2677 PetscErrorCode ierr; 2678 2679 PetscFunctionBegin; 2680 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2681 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2682 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2683 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2684 PetscFunctionReturn(0); 2685 } 2686 2687 /*@C 2688 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2689 2690 Collective 2691 2692 Input Parameters: 2693 + B - the matrix 2694 - data - the array (or NULL) 2695 2696 Notes: 2697 The data input variable is intended primarily for Fortran programmers 2698 who wish to allocate their own matrix memory space. Most users should 2699 need not call this routine. 2700 2701 Level: intermediate 2702 2703 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2704 2705 @*/ 2706 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2707 { 2708 PetscErrorCode ierr; 2709 2710 PetscFunctionBegin; 2711 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2712 PetscFunctionReturn(0); 2713 } 2714 2715 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2716 { 2717 Mat_SeqDense *b; 2718 PetscErrorCode ierr; 2719 2720 PetscFunctionBegin; 2721 B->preallocated = PETSC_TRUE; 2722 2723 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2724 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2725 2726 b = (Mat_SeqDense*)B->data; 2727 b->Mmax = B->rmap->n; 2728 b->Nmax = B->cmap->n; 2729 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2730 2731 ierr = PetscIntMultError(b->lda,b->Nmax,NULL);CHKERRQ(ierr); 2732 if (!data) { /* petsc-allocated storage */ 2733 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2734 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2735 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2736 2737 b->user_alloc = PETSC_FALSE; 2738 } else { /* user-allocated storage */ 2739 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2740 b->v = data; 2741 b->user_alloc = PETSC_TRUE; 2742 } 2743 B->assembled = PETSC_TRUE; 2744 PetscFunctionReturn(0); 2745 } 2746 2747 #if defined(PETSC_HAVE_ELEMENTAL) 2748 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2749 { 2750 Mat mat_elemental; 2751 PetscErrorCode ierr; 2752 const PetscScalar *array; 2753 PetscScalar *v_colwise; 2754 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2755 2756 PetscFunctionBegin; 2757 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2758 ierr = MatDenseGetArrayRead(A,&array);CHKERRQ(ierr); 2759 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2760 k = 0; 2761 for (j=0; j<N; j++) { 2762 cols[j] = j; 2763 for (i=0; i<M; i++) { 2764 v_colwise[j*M+i] = array[k++]; 2765 } 2766 } 2767 for (i=0; i<M; i++) { 2768 rows[i] = i; 2769 } 2770 ierr = MatDenseRestoreArrayRead(A,&array);CHKERRQ(ierr); 2771 2772 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2773 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2774 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2775 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2776 2777 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2778 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2779 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2780 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2781 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2782 2783 if (reuse == MAT_INPLACE_MATRIX) { 2784 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2785 } else { 2786 *newmat = mat_elemental; 2787 } 2788 PetscFunctionReturn(0); 2789 } 2790 #endif 2791 2792 /*@C 2793 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2794 2795 Input parameter: 2796 + A - the matrix 2797 - lda - the leading dimension 2798 2799 Notes: 2800 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2801 it asserts that the preallocation has a leading dimension (the LDA parameter 2802 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2803 2804 Level: intermediate 2805 2806 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2807 2808 @*/ 2809 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2810 { 2811 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2812 2813 PetscFunctionBegin; 2814 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); 2815 b->lda = lda; 2816 b->changelda = PETSC_FALSE; 2817 b->Mmax = PetscMax(b->Mmax,lda); 2818 PetscFunctionReturn(0); 2819 } 2820 2821 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2822 { 2823 PetscErrorCode ierr; 2824 PetscMPIInt size; 2825 2826 PetscFunctionBegin; 2827 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2828 if (size == 1) { 2829 if (scall == MAT_INITIAL_MATRIX) { 2830 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 2831 } else { 2832 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2833 } 2834 } else { 2835 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2836 } 2837 PetscFunctionReturn(0); 2838 } 2839 2840 PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v) 2841 { 2842 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2843 PetscErrorCode ierr; 2844 2845 PetscFunctionBegin; 2846 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec first"); 2847 if (!a->cvec) { 2848 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 2849 } 2850 a->vecinuse = col + 1; 2851 ierr = MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2852 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 2853 *v = a->cvec; 2854 PetscFunctionReturn(0); 2855 } 2856 2857 PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v) 2858 { 2859 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2860 PetscErrorCode ierr; 2861 2862 PetscFunctionBegin; 2863 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec first"); 2864 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 2865 a->vecinuse = 0; 2866 ierr = MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2867 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 2868 *v = NULL; 2869 PetscFunctionReturn(0); 2870 } 2871 2872 PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v) 2873 { 2874 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2875 PetscErrorCode ierr; 2876 2877 PetscFunctionBegin; 2878 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec first"); 2879 if (!a->cvec) { 2880 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 2881 } 2882 a->vecinuse = col + 1; 2883 ierr = MatDenseGetArrayRead(A,&a->ptrinuse);CHKERRQ(ierr); 2884 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 2885 ierr = VecLockReadPush(a->cvec);CHKERRQ(ierr); 2886 *v = a->cvec; 2887 PetscFunctionReturn(0); 2888 } 2889 2890 PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v) 2891 { 2892 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2893 PetscErrorCode ierr; 2894 2895 PetscFunctionBegin; 2896 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec first"); 2897 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 2898 a->vecinuse = 0; 2899 ierr = MatDenseRestoreArrayRead(A,&a->ptrinuse);CHKERRQ(ierr); 2900 ierr = VecLockReadPop(a->cvec);CHKERRQ(ierr); 2901 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 2902 *v = NULL; 2903 PetscFunctionReturn(0); 2904 } 2905 2906 PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v) 2907 { 2908 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2909 PetscErrorCode ierr; 2910 2911 PetscFunctionBegin; 2912 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec first"); 2913 if (!a->cvec) { 2914 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 2915 } 2916 a->vecinuse = col + 1; 2917 ierr = MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2918 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 2919 *v = a->cvec; 2920 PetscFunctionReturn(0); 2921 } 2922 2923 PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v) 2924 { 2925 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2926 PetscErrorCode ierr; 2927 2928 PetscFunctionBegin; 2929 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec first"); 2930 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 2931 a->vecinuse = 0; 2932 ierr = MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2933 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 2934 *v = NULL; 2935 PetscFunctionReturn(0); 2936 } 2937 2938 /*MC 2939 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2940 2941 Options Database Keys: 2942 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2943 2944 Level: beginner 2945 2946 .seealso: MatCreateSeqDense() 2947 2948 M*/ 2949 PetscErrorCode MatCreate_SeqDense(Mat B) 2950 { 2951 Mat_SeqDense *b; 2952 PetscErrorCode ierr; 2953 PetscMPIInt size; 2954 2955 PetscFunctionBegin; 2956 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2957 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2958 2959 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2960 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2961 B->data = (void*)b; 2962 2963 b->roworiented = PETSC_TRUE; 2964 2965 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);CHKERRQ(ierr); 2966 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2967 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2968 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 2969 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 2970 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2971 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2972 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2973 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2974 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2975 #if defined(PETSC_HAVE_ELEMENTAL) 2976 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 2977 #endif 2978 #if defined(PETSC_HAVE_CUDA) 2979 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);CHKERRQ(ierr); 2980 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 2981 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 2982 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 2983 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaijcusparse_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);CHKERRQ(ierr); 2984 #endif 2985 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2986 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);CHKERRQ(ierr); 2987 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 2988 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2989 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2990 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr); 2991 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqbaij_seqdense_C",MatMatMultSymbolic_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2992 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqbaij_seqdense_C",MatMatMultNumeric_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2993 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr); 2994 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqsbaij_seqdense_C",MatMatMultSymbolic_SeqSBAIJ_SeqDense);CHKERRQ(ierr); 2995 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqsbaij_seqdense_C",MatMatMultNumeric_SeqSBAIJ_SeqDense);CHKERRQ(ierr); 2996 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2997 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2998 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2999 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 3000 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 3001 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 3002 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_nest_seqdense_C",MatMatMultSymbolic_Nest_Dense);CHKERRQ(ierr); 3003 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_nest_seqdense_C",MatMatMultNumeric_Nest_Dense);CHKERRQ(ierr); 3004 3005 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 3006 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 3007 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 3008 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 3009 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 3010 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 3011 3012 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 3013 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 3014 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr); 3015 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr); 3016 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);CHKERRQ(ierr); 3017 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);CHKERRQ(ierr); 3018 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);CHKERRQ(ierr); 3019 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);CHKERRQ(ierr); 3020 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);CHKERRQ(ierr); 3021 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);CHKERRQ(ierr); 3022 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 3023 PetscFunctionReturn(0); 3024 } 3025 3026 /*@C 3027 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. 3028 3029 Not Collective 3030 3031 Input Parameter: 3032 + mat - a MATSEQDENSE or MATMPIDENSE matrix 3033 - col - column index 3034 3035 Output Parameter: 3036 . vals - pointer to the data 3037 3038 Level: intermediate 3039 3040 .seealso: MatDenseRestoreColumn() 3041 @*/ 3042 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 3043 { 3044 PetscErrorCode ierr; 3045 3046 PetscFunctionBegin; 3047 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr); 3048 PetscFunctionReturn(0); 3049 } 3050 3051 /*@C 3052 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 3053 3054 Not Collective 3055 3056 Input Parameter: 3057 . mat - a MATSEQDENSE or MATMPIDENSE matrix 3058 3059 Output Parameter: 3060 . vals - pointer to the data 3061 3062 Level: intermediate 3063 3064 .seealso: MatDenseGetColumn() 3065 @*/ 3066 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 3067 { 3068 PetscErrorCode ierr; 3069 3070 PetscFunctionBegin; 3071 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr); 3072 PetscFunctionReturn(0); 3073 } 3074 3075 /*@C 3076 MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec. 3077 3078 Collective 3079 3080 Input Parameter: 3081 + mat - the Mat object 3082 - col - the column index 3083 3084 Output Parameter: 3085 . v - the vector 3086 3087 Notes: 3088 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVec() when the vector is no longer needed. 3089 Use MatDenseGetColumnVecRead() to obtain read-only access or MatDenseGetColumnVecWrite() for write-only access. 3090 3091 Level: intermediate 3092 3093 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3094 @*/ 3095 PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v) 3096 { 3097 PetscErrorCode ierr; 3098 3099 PetscFunctionBegin; 3100 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3101 PetscValidType(A,1); 3102 PetscValidLogicalCollectiveInt(A,col,2); 3103 PetscValidPointer(v,3); 3104 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3105 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); 3106 ierr = PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3107 PetscFunctionReturn(0); 3108 } 3109 3110 /*@C 3111 MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec(). 3112 3113 Collective 3114 3115 Input Parameter: 3116 + mat - the Mat object 3117 . col - the column index 3118 - v - the Vec object 3119 3120 Level: intermediate 3121 3122 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3123 @*/ 3124 PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v) 3125 { 3126 PetscErrorCode ierr; 3127 3128 PetscFunctionBegin; 3129 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3130 PetscValidType(A,1); 3131 PetscValidLogicalCollectiveInt(A,col,2); 3132 PetscValidPointer(v,3); 3133 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3134 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); 3135 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3136 PetscFunctionReturn(0); 3137 } 3138 3139 /*@C 3140 MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec. 3141 3142 Collective 3143 3144 Input Parameter: 3145 + mat - the Mat object 3146 - col - the column index 3147 3148 Output Parameter: 3149 . v - the vector 3150 3151 Notes: 3152 The vector is owned by PETSc and users cannot modify it. 3153 Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed. 3154 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access. 3155 3156 Level: intermediate 3157 3158 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3159 @*/ 3160 PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v) 3161 { 3162 PetscErrorCode ierr; 3163 3164 PetscFunctionBegin; 3165 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3166 PetscValidType(A,1); 3167 PetscValidLogicalCollectiveInt(A,col,2); 3168 PetscValidPointer(v,3); 3169 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3170 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); 3171 ierr = PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3172 PetscFunctionReturn(0); 3173 } 3174 3175 /*@C 3176 MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead(). 3177 3178 Collective 3179 3180 Input Parameter: 3181 + mat - the Mat object 3182 . col - the column index 3183 - v - the Vec object 3184 3185 Level: intermediate 3186 3187 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite() 3188 @*/ 3189 PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v) 3190 { 3191 PetscErrorCode ierr; 3192 3193 PetscFunctionBegin; 3194 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3195 PetscValidType(A,1); 3196 PetscValidLogicalCollectiveInt(A,col,2); 3197 PetscValidPointer(v,3); 3198 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3199 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); 3200 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3201 PetscFunctionReturn(0); 3202 } 3203 3204 /*@C 3205 MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec. 3206 3207 Collective 3208 3209 Input Parameter: 3210 + mat - the Mat object 3211 - col - the column index 3212 3213 Output Parameter: 3214 . v - the vector 3215 3216 Notes: 3217 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed. 3218 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access. 3219 3220 Level: intermediate 3221 3222 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3223 @*/ 3224 PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v) 3225 { 3226 PetscErrorCode ierr; 3227 3228 PetscFunctionBegin; 3229 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3230 PetscValidType(A,1); 3231 PetscValidLogicalCollectiveInt(A,col,2); 3232 PetscValidPointer(v,3); 3233 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3234 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); 3235 ierr = PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3236 PetscFunctionReturn(0); 3237 } 3238 3239 /*@C 3240 MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite(). 3241 3242 Collective 3243 3244 Input Parameter: 3245 + mat - the Mat object 3246 . col - the column index 3247 - v - the Vec object 3248 3249 Level: intermediate 3250 3251 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead() 3252 @*/ 3253 PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v) 3254 { 3255 PetscErrorCode ierr; 3256 3257 PetscFunctionBegin; 3258 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3259 PetscValidType(A,1); 3260 PetscValidLogicalCollectiveInt(A,col,2); 3261 PetscValidPointer(v,3); 3262 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3263 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); 3264 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3265 PetscFunctionReturn(0); 3266 } 3267