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