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