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