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