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,"MatDenseGetArray_C",NULL);CHKERRQ(ierr); 1465 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);CHKERRQ(ierr); 1466 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);CHKERRQ(ierr); 1467 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr); 1468 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);CHKERRQ(ierr); 1469 #if defined(PETSC_HAVE_ELEMENTAL) 1470 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);CHKERRQ(ierr); 1471 #endif 1472 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1473 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1474 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1475 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1476 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1477 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1478 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1479 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1480 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);CHKERRQ(ierr); 1481 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);CHKERRQ(ierr); 1482 PetscFunctionReturn(0); 1483 } 1484 1485 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1486 { 1487 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1488 PetscErrorCode ierr; 1489 PetscInt k,j,m,n,M; 1490 PetscScalar *v,tmp; 1491 1492 PetscFunctionBegin; 1493 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1494 if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */ 1495 for (j=0; j<m; j++) { 1496 for (k=0; k<j; k++) { 1497 tmp = v[j + k*M]; 1498 v[j + k*M] = v[k + j*M]; 1499 v[k + j*M] = tmp; 1500 } 1501 } 1502 } else { /* out-of-place transpose */ 1503 Mat tmat; 1504 Mat_SeqDense *tmatd; 1505 PetscScalar *v2; 1506 PetscInt M2; 1507 1508 if (reuse != MAT_REUSE_MATRIX) { 1509 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1510 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1511 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1512 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1513 } else { 1514 tmat = *matout; 1515 } 1516 tmatd = (Mat_SeqDense*)tmat->data; 1517 v = mat->v; v2 = tmatd->v; M2 = tmatd->lda; 1518 for (j=0; j<n; j++) { 1519 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1520 } 1521 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1522 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1523 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat; 1524 else { 1525 ierr = MatHeaderMerge(A,&tmat);CHKERRQ(ierr); 1526 } 1527 } 1528 PetscFunctionReturn(0); 1529 } 1530 1531 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1532 { 1533 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1534 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1535 PetscInt i,j; 1536 PetscScalar *v1,*v2; 1537 1538 PetscFunctionBegin; 1539 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1540 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1541 for (i=0; i<A1->rmap->n; i++) { 1542 v1 = mat1->v+i; v2 = mat2->v+i; 1543 for (j=0; j<A1->cmap->n; j++) { 1544 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1545 v1 += mat1->lda; v2 += mat2->lda; 1546 } 1547 } 1548 *flg = PETSC_TRUE; 1549 PetscFunctionReturn(0); 1550 } 1551 1552 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1553 { 1554 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1555 PetscErrorCode ierr; 1556 PetscInt i,n,len; 1557 PetscScalar *x,zero = 0.0; 1558 1559 PetscFunctionBegin; 1560 ierr = VecSet(v,zero);CHKERRQ(ierr); 1561 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1562 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1563 len = PetscMin(A->rmap->n,A->cmap->n); 1564 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1565 for (i=0; i<len; i++) { 1566 x[i] = mat->v[i*mat->lda + i]; 1567 } 1568 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1569 PetscFunctionReturn(0); 1570 } 1571 1572 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1573 { 1574 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1575 const PetscScalar *l,*r; 1576 PetscScalar x,*v; 1577 PetscErrorCode ierr; 1578 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1579 1580 PetscFunctionBegin; 1581 if (ll) { 1582 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1583 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1584 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1585 for (i=0; i<m; i++) { 1586 x = l[i]; 1587 v = mat->v + i; 1588 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1589 } 1590 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1591 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1592 } 1593 if (rr) { 1594 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1595 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1596 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1597 for (i=0; i<n; i++) { 1598 x = r[i]; 1599 v = mat->v + i*mat->lda; 1600 for (j=0; j<m; j++) (*v++) *= x; 1601 } 1602 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1603 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1604 } 1605 PetscFunctionReturn(0); 1606 } 1607 1608 static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1609 { 1610 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1611 PetscScalar *v = mat->v; 1612 PetscReal sum = 0.0; 1613 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1614 PetscErrorCode ierr; 1615 1616 PetscFunctionBegin; 1617 if (type == NORM_FROBENIUS) { 1618 if (lda>m) { 1619 for (j=0; j<A->cmap->n; j++) { 1620 v = mat->v+j*lda; 1621 for (i=0; i<m; i++) { 1622 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1623 } 1624 } 1625 } else { 1626 #if defined(PETSC_USE_REAL___FP16) 1627 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1628 *nrm = BLASnrm2_(&cnt,v,&one); 1629 } 1630 #else 1631 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1632 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1633 } 1634 } 1635 *nrm = PetscSqrtReal(sum); 1636 #endif 1637 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1638 } else if (type == NORM_1) { 1639 *nrm = 0.0; 1640 for (j=0; j<A->cmap->n; j++) { 1641 v = mat->v + j*mat->lda; 1642 sum = 0.0; 1643 for (i=0; i<A->rmap->n; i++) { 1644 sum += PetscAbsScalar(*v); v++; 1645 } 1646 if (sum > *nrm) *nrm = sum; 1647 } 1648 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1649 } else if (type == NORM_INFINITY) { 1650 *nrm = 0.0; 1651 for (j=0; j<A->rmap->n; j++) { 1652 v = mat->v + j; 1653 sum = 0.0; 1654 for (i=0; i<A->cmap->n; i++) { 1655 sum += PetscAbsScalar(*v); v += mat->lda; 1656 } 1657 if (sum > *nrm) *nrm = sum; 1658 } 1659 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1660 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1661 PetscFunctionReturn(0); 1662 } 1663 1664 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1665 { 1666 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1667 PetscErrorCode ierr; 1668 1669 PetscFunctionBegin; 1670 switch (op) { 1671 case MAT_ROW_ORIENTED: 1672 aij->roworiented = flg; 1673 break; 1674 case MAT_NEW_NONZERO_LOCATIONS: 1675 case MAT_NEW_NONZERO_LOCATION_ERR: 1676 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1677 case MAT_NEW_DIAGONALS: 1678 case MAT_KEEP_NONZERO_PATTERN: 1679 case MAT_IGNORE_OFF_PROC_ENTRIES: 1680 case MAT_USE_HASH_TABLE: 1681 case MAT_IGNORE_ZERO_ENTRIES: 1682 case MAT_IGNORE_LOWER_TRIANGULAR: 1683 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1684 break; 1685 case MAT_SPD: 1686 case MAT_SYMMETRIC: 1687 case MAT_STRUCTURALLY_SYMMETRIC: 1688 case MAT_HERMITIAN: 1689 case MAT_SYMMETRY_ETERNAL: 1690 /* These options are handled directly by MatSetOption() */ 1691 break; 1692 default: 1693 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1694 } 1695 PetscFunctionReturn(0); 1696 } 1697 1698 static PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1699 { 1700 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1701 PetscErrorCode ierr; 1702 PetscInt lda=l->lda,m=A->rmap->n,j; 1703 1704 PetscFunctionBegin; 1705 if (lda>m) { 1706 for (j=0; j<A->cmap->n; j++) { 1707 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1708 } 1709 } else { 1710 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1711 } 1712 PetscFunctionReturn(0); 1713 } 1714 1715 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1716 { 1717 PetscErrorCode ierr; 1718 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1719 PetscInt m = l->lda, n = A->cmap->n, i,j; 1720 PetscScalar *slot,*bb; 1721 const PetscScalar *xx; 1722 1723 PetscFunctionBegin; 1724 #if defined(PETSC_USE_DEBUG) 1725 for (i=0; i<N; i++) { 1726 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1727 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); 1728 } 1729 #endif 1730 1731 /* fix right hand side if needed */ 1732 if (x && b) { 1733 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1734 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1735 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1736 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1737 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1738 } 1739 1740 for (i=0; i<N; i++) { 1741 slot = l->v + rows[i]; 1742 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1743 } 1744 if (diag != 0.0) { 1745 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1746 for (i=0; i<N; i++) { 1747 slot = l->v + (m+1)*rows[i]; 1748 *slot = diag; 1749 } 1750 } 1751 PetscFunctionReturn(0); 1752 } 1753 1754 static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1755 { 1756 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1757 1758 PetscFunctionBegin; 1759 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"); 1760 *array = mat->v; 1761 PetscFunctionReturn(0); 1762 } 1763 1764 static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1765 { 1766 PetscFunctionBegin; 1767 PetscFunctionReturn(0); 1768 } 1769 1770 /*@C 1771 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1772 1773 Logically Collective on Mat 1774 1775 Input Parameter: 1776 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1777 1778 Output Parameter: 1779 . array - pointer to the data 1780 1781 Level: intermediate 1782 1783 .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1784 @*/ 1785 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1786 { 1787 PetscErrorCode ierr; 1788 1789 PetscFunctionBegin; 1790 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1791 PetscFunctionReturn(0); 1792 } 1793 1794 /*@C 1795 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1796 1797 Logically Collective on Mat 1798 1799 Input Parameters: 1800 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1801 . array - pointer to the data 1802 1803 Level: intermediate 1804 1805 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1806 @*/ 1807 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1808 { 1809 PetscErrorCode ierr; 1810 1811 PetscFunctionBegin; 1812 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1813 if (array) *array = NULL; 1814 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1815 PetscFunctionReturn(0); 1816 } 1817 1818 /*@C 1819 MatDenseGetArrayRead - gives access to the array where the data for a SeqDense matrix is stored 1820 1821 Not Collective 1822 1823 Input Parameter: 1824 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1825 1826 Output Parameter: 1827 . array - pointer to the data 1828 1829 Level: intermediate 1830 1831 .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead() 1832 @*/ 1833 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 1834 { 1835 PetscErrorCode ierr; 1836 1837 PetscFunctionBegin; 1838 ierr = PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1839 PetscFunctionReturn(0); 1840 } 1841 1842 /*@C 1843 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1844 1845 Not Collective 1846 1847 Input Parameters: 1848 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1849 . array - pointer to the data 1850 1851 Level: intermediate 1852 1853 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray() 1854 @*/ 1855 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 1856 { 1857 PetscErrorCode ierr; 1858 1859 PetscFunctionBegin; 1860 ierr = PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1861 if (array) *array = NULL; 1862 PetscFunctionReturn(0); 1863 } 1864 1865 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1866 { 1867 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1868 PetscErrorCode ierr; 1869 PetscInt i,j,nrows,ncols; 1870 const PetscInt *irow,*icol; 1871 PetscScalar *av,*bv,*v = mat->v; 1872 Mat newmat; 1873 1874 PetscFunctionBegin; 1875 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1876 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1877 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1878 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1879 1880 /* Check submatrixcall */ 1881 if (scall == MAT_REUSE_MATRIX) { 1882 PetscInt n_cols,n_rows; 1883 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1884 if (n_rows != nrows || n_cols != ncols) { 1885 /* resize the result matrix to match number of requested rows/columns */ 1886 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1887 } 1888 newmat = *B; 1889 } else { 1890 /* Create and fill new matrix */ 1891 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1892 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1893 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1894 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1895 } 1896 1897 /* Now extract the data pointers and do the copy,column at a time */ 1898 bv = ((Mat_SeqDense*)newmat->data)->v; 1899 1900 for (i=0; i<ncols; i++) { 1901 av = v + mat->lda*icol[i]; 1902 for (j=0; j<nrows; j++) *bv++ = av[irow[j]]; 1903 } 1904 1905 /* Assemble the matrices so that the correct flags are set */ 1906 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1907 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1908 1909 /* Free work space */ 1910 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1911 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1912 *B = newmat; 1913 PetscFunctionReturn(0); 1914 } 1915 1916 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1917 { 1918 PetscErrorCode ierr; 1919 PetscInt i; 1920 1921 PetscFunctionBegin; 1922 if (scall == MAT_INITIAL_MATRIX) { 1923 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 1924 } 1925 1926 for (i=0; i<n; i++) { 1927 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1928 } 1929 PetscFunctionReturn(0); 1930 } 1931 1932 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1933 { 1934 PetscFunctionBegin; 1935 PetscFunctionReturn(0); 1936 } 1937 1938 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1939 { 1940 PetscFunctionBegin; 1941 PetscFunctionReturn(0); 1942 } 1943 1944 static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1945 { 1946 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 1947 PetscErrorCode ierr; 1948 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1949 1950 PetscFunctionBegin; 1951 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1952 if (A->ops->copy != B->ops->copy) { 1953 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1954 PetscFunctionReturn(0); 1955 } 1956 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1957 if (lda1>m || lda2>m) { 1958 for (j=0; j<n; j++) { 1959 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1960 } 1961 } else { 1962 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1963 } 1964 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 1965 PetscFunctionReturn(0); 1966 } 1967 1968 static PetscErrorCode MatSetUp_SeqDense(Mat A) 1969 { 1970 PetscErrorCode ierr; 1971 1972 PetscFunctionBegin; 1973 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1974 PetscFunctionReturn(0); 1975 } 1976 1977 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1978 { 1979 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1980 PetscInt i,nz = A->rmap->n*A->cmap->n; 1981 PetscScalar *aa = a->v; 1982 1983 PetscFunctionBegin; 1984 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1985 PetscFunctionReturn(0); 1986 } 1987 1988 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1989 { 1990 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1991 PetscInt i,nz = A->rmap->n*A->cmap->n; 1992 PetscScalar *aa = a->v; 1993 1994 PetscFunctionBegin; 1995 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1996 PetscFunctionReturn(0); 1997 } 1998 1999 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2000 { 2001 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2002 PetscInt i,nz = A->rmap->n*A->cmap->n; 2003 PetscScalar *aa = a->v; 2004 2005 PetscFunctionBegin; 2006 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 2007 PetscFunctionReturn(0); 2008 } 2009 2010 /* ----------------------------------------------------------------*/ 2011 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2012 { 2013 PetscErrorCode ierr; 2014 2015 PetscFunctionBegin; 2016 if (scall == MAT_INITIAL_MATRIX) { 2017 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2018 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2019 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2020 } 2021 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2022 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2023 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2024 PetscFunctionReturn(0); 2025 } 2026 2027 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2028 { 2029 PetscErrorCode ierr; 2030 PetscInt m=A->rmap->n,n=B->cmap->n; 2031 Mat Cmat; 2032 2033 PetscFunctionBegin; 2034 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); 2035 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2036 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2037 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 2038 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2039 2040 *C = Cmat; 2041 PetscFunctionReturn(0); 2042 } 2043 2044 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2045 { 2046 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2047 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2048 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2049 PetscBLASInt m,n,k; 2050 PetscScalar _DOne=1.0,_DZero=0.0; 2051 PetscErrorCode ierr; 2052 PetscBool flg; 2053 2054 PetscFunctionBegin; 2055 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);CHKERRQ(ierr); 2056 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Second matrix must be dense"); 2057 2058 /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/ 2059 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 2060 if (flg) { 2061 C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense; 2062 ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr); 2063 PetscFunctionReturn(0); 2064 } 2065 2066 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);CHKERRQ(ierr); 2067 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense"); 2068 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2069 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2070 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2071 if (!m || !n || !k) PetscFunctionReturn(0); 2072 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2073 PetscFunctionReturn(0); 2074 } 2075 2076 PetscErrorCode MatMatTransposeMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2077 { 2078 PetscErrorCode ierr; 2079 2080 PetscFunctionBegin; 2081 if (scall == MAT_INITIAL_MATRIX) { 2082 ierr = PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2083 ierr = MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2084 ierr = PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2085 } 2086 ierr = PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 2087 ierr = MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2088 ierr = PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 2089 PetscFunctionReturn(0); 2090 } 2091 2092 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2093 { 2094 PetscErrorCode ierr; 2095 PetscInt m=A->rmap->n,n=B->rmap->n; 2096 Mat Cmat; 2097 2098 PetscFunctionBegin; 2099 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); 2100 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2101 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2102 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 2103 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2104 2105 Cmat->assembled = PETSC_TRUE; 2106 2107 *C = Cmat; 2108 PetscFunctionReturn(0); 2109 } 2110 2111 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2112 { 2113 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2114 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2115 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2116 PetscBLASInt m,n,k; 2117 PetscScalar _DOne=1.0,_DZero=0.0; 2118 PetscErrorCode ierr; 2119 2120 PetscFunctionBegin; 2121 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2122 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2123 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2124 if (!m || !n || !k) PetscFunctionReturn(0); 2125 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2126 PetscFunctionReturn(0); 2127 } 2128 2129 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2130 { 2131 PetscErrorCode ierr; 2132 2133 PetscFunctionBegin; 2134 if (scall == MAT_INITIAL_MATRIX) { 2135 ierr = PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2136 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2137 ierr = PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2138 } 2139 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2140 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2141 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2142 PetscFunctionReturn(0); 2143 } 2144 2145 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2146 { 2147 PetscErrorCode ierr; 2148 PetscInt m=A->cmap->n,n=B->cmap->n; 2149 Mat Cmat; 2150 2151 PetscFunctionBegin; 2152 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); 2153 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2154 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2155 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 2156 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2157 2158 Cmat->assembled = PETSC_TRUE; 2159 2160 *C = Cmat; 2161 PetscFunctionReturn(0); 2162 } 2163 2164 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2165 { 2166 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2167 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2168 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2169 PetscBLASInt m,n,k; 2170 PetscScalar _DOne=1.0,_DZero=0.0; 2171 PetscErrorCode ierr; 2172 2173 PetscFunctionBegin; 2174 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2175 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2176 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 2177 if (!m || !n || !k) PetscFunctionReturn(0); 2178 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2179 PetscFunctionReturn(0); 2180 } 2181 2182 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2183 { 2184 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2185 PetscErrorCode ierr; 2186 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2187 PetscScalar *x; 2188 MatScalar *aa = a->v; 2189 2190 PetscFunctionBegin; 2191 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2192 2193 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2194 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2195 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2196 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2197 for (i=0; i<m; i++) { 2198 x[i] = aa[i]; if (idx) idx[i] = 0; 2199 for (j=1; j<n; j++) { 2200 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 2201 } 2202 } 2203 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2204 PetscFunctionReturn(0); 2205 } 2206 2207 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2208 { 2209 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2210 PetscErrorCode ierr; 2211 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2212 PetscScalar *x; 2213 PetscReal atmp; 2214 MatScalar *aa = a->v; 2215 2216 PetscFunctionBegin; 2217 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2218 2219 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2220 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2221 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2222 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2223 for (i=0; i<m; i++) { 2224 x[i] = PetscAbsScalar(aa[i]); 2225 for (j=1; j<n; j++) { 2226 atmp = PetscAbsScalar(aa[i+m*j]); 2227 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2228 } 2229 } 2230 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2231 PetscFunctionReturn(0); 2232 } 2233 2234 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2235 { 2236 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2237 PetscErrorCode ierr; 2238 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2239 PetscScalar *x; 2240 MatScalar *aa = a->v; 2241 2242 PetscFunctionBegin; 2243 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2244 2245 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2246 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2247 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2248 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2249 for (i=0; i<m; i++) { 2250 x[i] = aa[i]; if (idx) idx[i] = 0; 2251 for (j=1; j<n; j++) { 2252 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 2253 } 2254 } 2255 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2256 PetscFunctionReturn(0); 2257 } 2258 2259 static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2260 { 2261 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2262 PetscErrorCode ierr; 2263 PetscScalar *x; 2264 2265 PetscFunctionBegin; 2266 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2267 2268 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2269 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 2270 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2271 PetscFunctionReturn(0); 2272 } 2273 2274 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2275 { 2276 PetscErrorCode ierr; 2277 PetscInt i,j,m,n; 2278 PetscScalar *a; 2279 2280 PetscFunctionBegin; 2281 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2282 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 2283 ierr = MatDenseGetArray(A,&a);CHKERRQ(ierr); 2284 if (type == NORM_2) { 2285 for (i=0; i<n; i++) { 2286 for (j=0; j<m; j++) { 2287 norms[i] += PetscAbsScalar(a[j]*a[j]); 2288 } 2289 a += m; 2290 } 2291 } else if (type == NORM_1) { 2292 for (i=0; i<n; i++) { 2293 for (j=0; j<m; j++) { 2294 norms[i] += PetscAbsScalar(a[j]); 2295 } 2296 a += m; 2297 } 2298 } else if (type == NORM_INFINITY) { 2299 for (i=0; i<n; i++) { 2300 for (j=0; j<m; j++) { 2301 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2302 } 2303 a += m; 2304 } 2305 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2306 ierr = MatDenseRestoreArray(A,&a);CHKERRQ(ierr); 2307 if (type == NORM_2) { 2308 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2309 } 2310 PetscFunctionReturn(0); 2311 } 2312 2313 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2314 { 2315 PetscErrorCode ierr; 2316 PetscScalar *a; 2317 PetscInt m,n,i; 2318 2319 PetscFunctionBegin; 2320 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2321 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2322 for (i=0; i<m*n; i++) { 2323 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 2324 } 2325 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2326 PetscFunctionReturn(0); 2327 } 2328 2329 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2330 { 2331 PetscFunctionBegin; 2332 *missing = PETSC_FALSE; 2333 PetscFunctionReturn(0); 2334 } 2335 2336 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2337 { 2338 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2339 2340 PetscFunctionBegin; 2341 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2342 *vals = a->v+col*a->lda; 2343 PetscFunctionReturn(0); 2344 } 2345 2346 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2347 { 2348 PetscFunctionBegin; 2349 *vals = 0; /* user cannot accidently use the array later */ 2350 PetscFunctionReturn(0); 2351 } 2352 2353 /* -------------------------------------------------------------------*/ 2354 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2355 MatGetRow_SeqDense, 2356 MatRestoreRow_SeqDense, 2357 MatMult_SeqDense, 2358 /* 4*/ MatMultAdd_SeqDense, 2359 MatMultTranspose_SeqDense, 2360 MatMultTransposeAdd_SeqDense, 2361 0, 2362 0, 2363 0, 2364 /* 10*/ 0, 2365 MatLUFactor_SeqDense, 2366 MatCholeskyFactor_SeqDense, 2367 MatSOR_SeqDense, 2368 MatTranspose_SeqDense, 2369 /* 15*/ MatGetInfo_SeqDense, 2370 MatEqual_SeqDense, 2371 MatGetDiagonal_SeqDense, 2372 MatDiagonalScale_SeqDense, 2373 MatNorm_SeqDense, 2374 /* 20*/ MatAssemblyBegin_SeqDense, 2375 MatAssemblyEnd_SeqDense, 2376 MatSetOption_SeqDense, 2377 MatZeroEntries_SeqDense, 2378 /* 24*/ MatZeroRows_SeqDense, 2379 0, 2380 0, 2381 0, 2382 0, 2383 /* 29*/ MatSetUp_SeqDense, 2384 0, 2385 0, 2386 0, 2387 0, 2388 /* 34*/ MatDuplicate_SeqDense, 2389 0, 2390 0, 2391 0, 2392 0, 2393 /* 39*/ MatAXPY_SeqDense, 2394 MatCreateSubMatrices_SeqDense, 2395 0, 2396 MatGetValues_SeqDense, 2397 MatCopy_SeqDense, 2398 /* 44*/ MatGetRowMax_SeqDense, 2399 MatScale_SeqDense, 2400 MatShift_Basic, 2401 0, 2402 MatZeroRowsColumns_SeqDense, 2403 /* 49*/ MatSetRandom_SeqDense, 2404 0, 2405 0, 2406 0, 2407 0, 2408 /* 54*/ 0, 2409 0, 2410 0, 2411 0, 2412 0, 2413 /* 59*/ 0, 2414 MatDestroy_SeqDense, 2415 MatView_SeqDense, 2416 0, 2417 0, 2418 /* 64*/ 0, 2419 0, 2420 0, 2421 0, 2422 0, 2423 /* 69*/ MatGetRowMaxAbs_SeqDense, 2424 0, 2425 0, 2426 0, 2427 0, 2428 /* 74*/ 0, 2429 0, 2430 0, 2431 0, 2432 0, 2433 /* 79*/ 0, 2434 0, 2435 0, 2436 0, 2437 /* 83*/ MatLoad_SeqDense, 2438 0, 2439 MatIsHermitian_SeqDense, 2440 0, 2441 0, 2442 0, 2443 /* 89*/ MatMatMult_SeqDense_SeqDense, 2444 MatMatMultSymbolic_SeqDense_SeqDense, 2445 MatMatMultNumeric_SeqDense_SeqDense, 2446 MatPtAP_SeqDense_SeqDense, 2447 MatPtAPSymbolic_SeqDense_SeqDense, 2448 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense, 2449 MatMatTransposeMult_SeqDense_SeqDense, 2450 MatMatTransposeMultSymbolic_SeqDense_SeqDense, 2451 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2452 0, 2453 /* 99*/ 0, 2454 0, 2455 0, 2456 MatConjugate_SeqDense, 2457 0, 2458 /*104*/ 0, 2459 MatRealPart_SeqDense, 2460 MatImaginaryPart_SeqDense, 2461 0, 2462 0, 2463 /*109*/ 0, 2464 0, 2465 MatGetRowMin_SeqDense, 2466 MatGetColumnVector_SeqDense, 2467 MatMissingDiagonal_SeqDense, 2468 /*114*/ 0, 2469 0, 2470 0, 2471 0, 2472 0, 2473 /*119*/ 0, 2474 0, 2475 0, 2476 0, 2477 0, 2478 /*124*/ 0, 2479 MatGetColumnNorms_SeqDense, 2480 0, 2481 0, 2482 0, 2483 /*129*/ 0, 2484 MatTransposeMatMult_SeqDense_SeqDense, 2485 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2486 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2487 0, 2488 /*134*/ 0, 2489 0, 2490 0, 2491 0, 2492 0, 2493 /*139*/ 0, 2494 0, 2495 0, 2496 0, 2497 0, 2498 /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqDense 2499 }; 2500 2501 /*@C 2502 MatCreateSeqDense - Creates a sequential dense matrix that 2503 is stored in column major order (the usual Fortran 77 manner). Many 2504 of the matrix operations use the BLAS and LAPACK routines. 2505 2506 Collective on MPI_Comm 2507 2508 Input Parameters: 2509 + comm - MPI communicator, set to PETSC_COMM_SELF 2510 . m - number of rows 2511 . n - number of columns 2512 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2513 to control all matrix memory allocation. 2514 2515 Output Parameter: 2516 . A - the matrix 2517 2518 Notes: 2519 The data input variable is intended primarily for Fortran programmers 2520 who wish to allocate their own matrix memory space. Most users should 2521 set data=NULL. 2522 2523 Level: intermediate 2524 2525 .keywords: dense, matrix, LAPACK, BLAS 2526 2527 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2528 @*/ 2529 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2530 { 2531 PetscErrorCode ierr; 2532 2533 PetscFunctionBegin; 2534 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2535 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2536 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2537 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2538 PetscFunctionReturn(0); 2539 } 2540 2541 /*@C 2542 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2543 2544 Collective on MPI_Comm 2545 2546 Input Parameters: 2547 + B - the matrix 2548 - data - the array (or NULL) 2549 2550 Notes: 2551 The data input variable is intended primarily for Fortran programmers 2552 who wish to allocate their own matrix memory space. Most users should 2553 need not call this routine. 2554 2555 Level: intermediate 2556 2557 .keywords: dense, matrix, LAPACK, BLAS 2558 2559 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2560 2561 @*/ 2562 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2563 { 2564 PetscErrorCode ierr; 2565 2566 PetscFunctionBegin; 2567 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2568 PetscFunctionReturn(0); 2569 } 2570 2571 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2572 { 2573 Mat_SeqDense *b; 2574 PetscErrorCode ierr; 2575 2576 PetscFunctionBegin; 2577 B->preallocated = PETSC_TRUE; 2578 2579 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2580 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2581 2582 b = (Mat_SeqDense*)B->data; 2583 b->Mmax = B->rmap->n; 2584 b->Nmax = B->cmap->n; 2585 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2586 2587 ierr = PetscIntMultError(b->lda,b->Nmax,NULL);CHKERRQ(ierr); 2588 if (!data) { /* petsc-allocated storage */ 2589 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2590 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2591 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2592 2593 b->user_alloc = PETSC_FALSE; 2594 } else { /* user-allocated storage */ 2595 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2596 b->v = data; 2597 b->user_alloc = PETSC_TRUE; 2598 } 2599 B->assembled = PETSC_TRUE; 2600 PetscFunctionReturn(0); 2601 } 2602 2603 #if defined(PETSC_HAVE_ELEMENTAL) 2604 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2605 { 2606 Mat mat_elemental; 2607 PetscErrorCode ierr; 2608 PetscScalar *array,*v_colwise; 2609 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2610 2611 PetscFunctionBegin; 2612 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2613 ierr = MatDenseGetArray(A,&array);CHKERRQ(ierr); 2614 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2615 k = 0; 2616 for (j=0; j<N; j++) { 2617 cols[j] = j; 2618 for (i=0; i<M; i++) { 2619 v_colwise[j*M+i] = array[k++]; 2620 } 2621 } 2622 for (i=0; i<M; i++) { 2623 rows[i] = i; 2624 } 2625 ierr = MatDenseRestoreArray(A,&array);CHKERRQ(ierr); 2626 2627 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2628 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2629 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2630 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2631 2632 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2633 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2634 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2635 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2636 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2637 2638 if (reuse == MAT_INPLACE_MATRIX) { 2639 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2640 } else { 2641 *newmat = mat_elemental; 2642 } 2643 PetscFunctionReturn(0); 2644 } 2645 #endif 2646 2647 /*@C 2648 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2649 2650 Input parameter: 2651 + A - the matrix 2652 - lda - the leading dimension 2653 2654 Notes: 2655 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2656 it asserts that the preallocation has a leading dimension (the LDA parameter 2657 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2658 2659 Level: intermediate 2660 2661 .keywords: dense, matrix, LAPACK, BLAS 2662 2663 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2664 2665 @*/ 2666 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2667 { 2668 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2669 2670 PetscFunctionBegin; 2671 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); 2672 b->lda = lda; 2673 b->changelda = PETSC_FALSE; 2674 b->Mmax = PetscMax(b->Mmax,lda); 2675 PetscFunctionReturn(0); 2676 } 2677 2678 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2679 { 2680 PetscErrorCode ierr; 2681 PetscMPIInt size; 2682 2683 PetscFunctionBegin; 2684 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2685 if (size == 1) { 2686 if (scall == MAT_INITIAL_MATRIX) { 2687 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 2688 } else { 2689 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2690 } 2691 } else { 2692 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2693 } 2694 PetscFunctionReturn(0); 2695 } 2696 2697 /*MC 2698 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2699 2700 Options Database Keys: 2701 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2702 2703 Level: beginner 2704 2705 .seealso: MatCreateSeqDense() 2706 2707 M*/ 2708 2709 PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B) 2710 { 2711 Mat_SeqDense *b; 2712 PetscErrorCode ierr; 2713 PetscMPIInt size; 2714 2715 PetscFunctionBegin; 2716 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2717 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2718 2719 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2720 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2721 B->data = (void*)b; 2722 2723 b->roworiented = PETSC_TRUE; 2724 2725 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2726 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2727 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 2728 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 2729 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2730 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2731 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2732 #if defined(PETSC_HAVE_ELEMENTAL) 2733 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 2734 #endif 2735 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2736 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2737 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2738 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2739 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2740 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2741 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2742 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2743 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2744 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijsell_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2745 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2746 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2747 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijsell_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2748 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2749 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2750 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2751 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2752 2753 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2754 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2755 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2756 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2757 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2758 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2759 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijsell_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2760 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2761 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2762 2763 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2764 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2765 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2766 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr); 2767 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr); 2768 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2769 PetscFunctionReturn(0); 2770 } 2771 2772 /*@C 2773 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. 2774 2775 Not Collective 2776 2777 Input Parameter: 2778 + mat - a MATSEQDENSE or MATMPIDENSE matrix 2779 - col - column index 2780 2781 Output Parameter: 2782 . vals - pointer to the data 2783 2784 Level: intermediate 2785 2786 .seealso: MatDenseRestoreColumn() 2787 @*/ 2788 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 2789 { 2790 PetscErrorCode ierr; 2791 2792 PetscFunctionBegin; 2793 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr); 2794 PetscFunctionReturn(0); 2795 } 2796 2797 /*@C 2798 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 2799 2800 Not Collective 2801 2802 Input Parameter: 2803 . mat - a MATSEQDENSE or MATMPIDENSE matrix 2804 2805 Output Parameter: 2806 . vals - pointer to the data 2807 2808 Level: intermediate 2809 2810 .seealso: MatDenseGetColumn() 2811 @*/ 2812 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 2813 { 2814 PetscErrorCode ierr; 2815 2816 PetscFunctionBegin; 2817 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr); 2818 PetscFunctionReturn(0); 2819 } 2820