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