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 PetscBool isbinary; 1085 1086 PetscFunctionBegin; 1087 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1088 if (!isbinary) SETERRQ2(PetscObjectComm((PetscObject)newmat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newmat)->type_name); 1089 1090 /* force binary viewer to load .info file if it has not yet done so */ 1091 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1092 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 1093 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1094 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 1095 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1096 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1097 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 1098 M = header[1]; N = header[2]; nz = header[3]; 1099 1100 /* set global size if not set already*/ 1101 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 1102 ierr = MatSetSizes(newmat,M,N,M,N);CHKERRQ(ierr); 1103 } else { 1104 /* if sizes and type are already set, check if the vector global sizes are correct */ 1105 ierr = MatGetSize(newmat,&grows,&gcols);CHKERRQ(ierr); 1106 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); 1107 } 1108 a = (Mat_SeqDense*)newmat->data; 1109 if (!a->user_alloc) { 1110 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1111 } 1112 1113 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 1114 a = (Mat_SeqDense*)newmat->data; 1115 v = a->v; 1116 /* Allocate some temp space to read in the values and then flip them 1117 from row major to column major */ 1118 ierr = PetscMalloc1(M*N > 0 ? M*N : 1,&w);CHKERRQ(ierr); 1119 /* read in nonzero values */ 1120 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 1121 /* now flip the values and store them in the matrix*/ 1122 for (j=0; j<N; j++) { 1123 for (i=0; i<M; i++) { 1124 *v++ =w[i*N+j]; 1125 } 1126 } 1127 ierr = PetscFree(w);CHKERRQ(ierr); 1128 } else { 1129 /* read row lengths */ 1130 ierr = PetscMalloc1(M+1,&rowlengths);CHKERRQ(ierr); 1131 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1132 1133 a = (Mat_SeqDense*)newmat->data; 1134 v = a->v; 1135 1136 /* read column indices and nonzeros */ 1137 ierr = PetscMalloc1(nz+1,&scols);CHKERRQ(ierr); 1138 cols = scols; 1139 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 1140 ierr = PetscMalloc1(nz+1,&svals);CHKERRQ(ierr); 1141 vals = svals; 1142 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 1143 1144 /* insert into matrix */ 1145 for (i=0; i<M; i++) { 1146 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 1147 svals += rowlengths[i]; scols += rowlengths[i]; 1148 } 1149 ierr = PetscFree(vals);CHKERRQ(ierr); 1150 ierr = PetscFree(cols);CHKERRQ(ierr); 1151 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1152 } 1153 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1154 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1155 PetscFunctionReturn(0); 1156 } 1157 1158 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 1159 { 1160 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1161 PetscErrorCode ierr; 1162 PetscInt i,j; 1163 const char *name; 1164 PetscScalar *v; 1165 PetscViewerFormat format; 1166 #if defined(PETSC_USE_COMPLEX) 1167 PetscBool allreal = PETSC_TRUE; 1168 #endif 1169 1170 PetscFunctionBegin; 1171 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1172 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 1173 PetscFunctionReturn(0); /* do nothing for now */ 1174 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 1175 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1176 for (i=0; i<A->rmap->n; i++) { 1177 v = a->v + i; 1178 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 1179 for (j=0; j<A->cmap->n; j++) { 1180 #if defined(PETSC_USE_COMPLEX) 1181 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 1182 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1183 } else if (PetscRealPart(*v)) { 1184 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));CHKERRQ(ierr); 1185 } 1186 #else 1187 if (*v) { 1188 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);CHKERRQ(ierr); 1189 } 1190 #endif 1191 v += a->lda; 1192 } 1193 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1194 } 1195 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1196 } else { 1197 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1198 #if defined(PETSC_USE_COMPLEX) 1199 /* determine if matrix has all real values */ 1200 v = a->v; 1201 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 1202 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;} 1203 } 1204 #endif 1205 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1206 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 1207 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1208 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1209 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 1210 } 1211 1212 for (i=0; i<A->rmap->n; i++) { 1213 v = a->v + i; 1214 for (j=0; j<A->cmap->n; j++) { 1215 #if defined(PETSC_USE_COMPLEX) 1216 if (allreal) { 1217 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));CHKERRQ(ierr); 1218 } else { 1219 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1220 } 1221 #else 1222 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);CHKERRQ(ierr); 1223 #endif 1224 v += a->lda; 1225 } 1226 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1227 } 1228 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1229 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 1230 } 1231 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1232 } 1233 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1234 PetscFunctionReturn(0); 1235 } 1236 1237 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 1238 { 1239 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1240 PetscErrorCode ierr; 1241 int fd; 1242 PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n; 1243 PetscScalar *v,*anonz,*vals; 1244 PetscViewerFormat format; 1245 1246 PetscFunctionBegin; 1247 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1248 1249 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1250 if (format == PETSC_VIEWER_NATIVE) { 1251 /* store the matrix as a dense matrix */ 1252 ierr = PetscMalloc1(4,&col_lens);CHKERRQ(ierr); 1253 1254 col_lens[0] = MAT_FILE_CLASSID; 1255 col_lens[1] = m; 1256 col_lens[2] = n; 1257 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 1258 1259 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1260 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1261 1262 /* write out matrix, by rows */ 1263 ierr = PetscMalloc1(m*n+1,&vals);CHKERRQ(ierr); 1264 v = a->v; 1265 for (j=0; j<n; j++) { 1266 for (i=0; i<m; i++) { 1267 vals[j + i*n] = *v++; 1268 } 1269 } 1270 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1271 ierr = PetscFree(vals);CHKERRQ(ierr); 1272 } else { 1273 ierr = PetscMalloc1(4+nz,&col_lens);CHKERRQ(ierr); 1274 1275 col_lens[0] = MAT_FILE_CLASSID; 1276 col_lens[1] = m; 1277 col_lens[2] = n; 1278 col_lens[3] = nz; 1279 1280 /* store lengths of each row and write (including header) to file */ 1281 for (i=0; i<m; i++) col_lens[4+i] = n; 1282 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1283 1284 /* Possibly should write in smaller increments, not whole matrix at once? */ 1285 /* store column indices (zero start index) */ 1286 ict = 0; 1287 for (i=0; i<m; i++) { 1288 for (j=0; j<n; j++) col_lens[ict++] = j; 1289 } 1290 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 1291 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1292 1293 /* store nonzero values */ 1294 ierr = PetscMalloc1(nz+1,&anonz);CHKERRQ(ierr); 1295 ict = 0; 1296 for (i=0; i<m; i++) { 1297 v = a->v + i; 1298 for (j=0; j<n; j++) { 1299 anonz[ict++] = *v; v += a->lda; 1300 } 1301 } 1302 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1303 ierr = PetscFree(anonz);CHKERRQ(ierr); 1304 } 1305 PetscFunctionReturn(0); 1306 } 1307 1308 #include <petscdraw.h> 1309 static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1310 { 1311 Mat A = (Mat) Aa; 1312 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1313 PetscErrorCode ierr; 1314 PetscInt m = A->rmap->n,n = A->cmap->n,i,j; 1315 int color = PETSC_DRAW_WHITE; 1316 PetscScalar *v = a->v; 1317 PetscViewer viewer; 1318 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 1319 PetscViewerFormat format; 1320 1321 PetscFunctionBegin; 1322 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1323 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1324 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1325 1326 /* Loop over matrix elements drawing boxes */ 1327 1328 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1329 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1330 /* Blue for negative and Red for positive */ 1331 for (j = 0; j < n; j++) { 1332 x_l = j; x_r = x_l + 1.0; 1333 for (i = 0; i < m; i++) { 1334 y_l = m - i - 1.0; 1335 y_r = y_l + 1.0; 1336 if (PetscRealPart(v[j*m+i]) > 0.) { 1337 color = PETSC_DRAW_RED; 1338 } else if (PetscRealPart(v[j*m+i]) < 0.) { 1339 color = PETSC_DRAW_BLUE; 1340 } else { 1341 continue; 1342 } 1343 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1344 } 1345 } 1346 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1347 } else { 1348 /* use contour shading to indicate magnitude of values */ 1349 /* first determine max of all nonzero values */ 1350 PetscReal minv = 0.0, maxv = 0.0; 1351 PetscDraw popup; 1352 1353 for (i=0; i < m*n; i++) { 1354 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1355 } 1356 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1357 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1358 ierr = PetscDrawScalePopup(popup,minv,maxv);CHKERRQ(ierr); 1359 1360 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1361 for (j=0; j<n; j++) { 1362 x_l = j; 1363 x_r = x_l + 1.0; 1364 for (i=0; i<m; i++) { 1365 y_l = m - i - 1.0; 1366 y_r = y_l + 1.0; 1367 color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv); 1368 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1369 } 1370 } 1371 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1372 } 1373 PetscFunctionReturn(0); 1374 } 1375 1376 static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1377 { 1378 PetscDraw draw; 1379 PetscBool isnull; 1380 PetscReal xr,yr,xl,yl,h,w; 1381 PetscErrorCode ierr; 1382 1383 PetscFunctionBegin; 1384 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1385 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1386 if (isnull) PetscFunctionReturn(0); 1387 1388 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1389 xr += w; yr += h; xl = -w; yl = -h; 1390 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1391 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1392 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1393 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 1394 ierr = PetscDrawSave(draw);CHKERRQ(ierr); 1395 PetscFunctionReturn(0); 1396 } 1397 1398 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1399 { 1400 PetscErrorCode ierr; 1401 PetscBool iascii,isbinary,isdraw; 1402 1403 PetscFunctionBegin; 1404 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1405 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1406 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1407 1408 if (iascii) { 1409 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1410 } else if (isbinary) { 1411 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1412 } else if (isdraw) { 1413 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1414 } 1415 PetscFunctionReturn(0); 1416 } 1417 1418 static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar array[]) 1419 { 1420 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1421 1422 PetscFunctionBegin; 1423 a->unplacedarray = a->v; 1424 a->unplaced_user_alloc = a->user_alloc; 1425 a->v = (PetscScalar*) array; 1426 PetscFunctionReturn(0); 1427 } 1428 1429 static PetscErrorCode MatDenseResetArray_SeqDense(Mat A) 1430 { 1431 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1432 1433 PetscFunctionBegin; 1434 a->v = a->unplacedarray; 1435 a->user_alloc = a->unplaced_user_alloc; 1436 a->unplacedarray = NULL; 1437 PetscFunctionReturn(0); 1438 } 1439 1440 static PetscErrorCode MatDestroy_SeqDense(Mat mat) 1441 { 1442 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1443 PetscErrorCode ierr; 1444 1445 PetscFunctionBegin; 1446 #if defined(PETSC_USE_LOG) 1447 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1448 #endif 1449 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1450 ierr = PetscFree(l->fwork);CHKERRQ(ierr); 1451 ierr = MatDestroy(&l->ptapwork);CHKERRQ(ierr); 1452 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1453 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1454 1455 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1456 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);CHKERRQ(ierr); 1457 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);CHKERRQ(ierr); 1458 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);CHKERRQ(ierr); 1459 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr); 1460 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);CHKERRQ(ierr); 1461 #if defined(PETSC_HAVE_ELEMENTAL) 1462 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);CHKERRQ(ierr); 1463 #endif 1464 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1465 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1466 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1467 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1468 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1469 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1470 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1471 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1472 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);CHKERRQ(ierr); 1473 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);CHKERRQ(ierr); 1474 PetscFunctionReturn(0); 1475 } 1476 1477 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1478 { 1479 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1480 PetscErrorCode ierr; 1481 PetscInt k,j,m,n,M; 1482 PetscScalar *v,tmp; 1483 1484 PetscFunctionBegin; 1485 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1486 if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */ 1487 for (j=0; j<m; j++) { 1488 for (k=0; k<j; k++) { 1489 tmp = v[j + k*M]; 1490 v[j + k*M] = v[k + j*M]; 1491 v[k + j*M] = tmp; 1492 } 1493 } 1494 } else { /* out-of-place transpose */ 1495 Mat tmat; 1496 Mat_SeqDense *tmatd; 1497 PetscScalar *v2; 1498 PetscInt M2; 1499 1500 if (reuse != MAT_REUSE_MATRIX) { 1501 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1502 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1503 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1504 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1505 } else { 1506 tmat = *matout; 1507 } 1508 tmatd = (Mat_SeqDense*)tmat->data; 1509 v = mat->v; v2 = tmatd->v; M2 = tmatd->lda; 1510 for (j=0; j<n; j++) { 1511 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1512 } 1513 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1514 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1515 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat; 1516 else { 1517 ierr = MatHeaderMerge(A,&tmat);CHKERRQ(ierr); 1518 } 1519 } 1520 PetscFunctionReturn(0); 1521 } 1522 1523 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1524 { 1525 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1526 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1527 PetscInt i,j; 1528 PetscScalar *v1,*v2; 1529 1530 PetscFunctionBegin; 1531 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1532 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1533 for (i=0; i<A1->rmap->n; i++) { 1534 v1 = mat1->v+i; v2 = mat2->v+i; 1535 for (j=0; j<A1->cmap->n; j++) { 1536 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1537 v1 += mat1->lda; v2 += mat2->lda; 1538 } 1539 } 1540 *flg = PETSC_TRUE; 1541 PetscFunctionReturn(0); 1542 } 1543 1544 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1545 { 1546 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1547 PetscErrorCode ierr; 1548 PetscInt i,n,len; 1549 PetscScalar *x,zero = 0.0; 1550 1551 PetscFunctionBegin; 1552 ierr = VecSet(v,zero);CHKERRQ(ierr); 1553 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1554 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1555 len = PetscMin(A->rmap->n,A->cmap->n); 1556 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1557 for (i=0; i<len; i++) { 1558 x[i] = mat->v[i*mat->lda + i]; 1559 } 1560 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1561 PetscFunctionReturn(0); 1562 } 1563 1564 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1565 { 1566 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1567 const PetscScalar *l,*r; 1568 PetscScalar x,*v; 1569 PetscErrorCode ierr; 1570 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1571 1572 PetscFunctionBegin; 1573 if (ll) { 1574 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1575 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1576 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1577 for (i=0; i<m; i++) { 1578 x = l[i]; 1579 v = mat->v + i; 1580 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1581 } 1582 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1583 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1584 } 1585 if (rr) { 1586 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1587 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1588 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1589 for (i=0; i<n; i++) { 1590 x = r[i]; 1591 v = mat->v + i*mat->lda; 1592 for (j=0; j<m; j++) (*v++) *= x; 1593 } 1594 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1595 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1596 } 1597 PetscFunctionReturn(0); 1598 } 1599 1600 static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1601 { 1602 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1603 PetscScalar *v = mat->v; 1604 PetscReal sum = 0.0; 1605 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1606 PetscErrorCode ierr; 1607 1608 PetscFunctionBegin; 1609 if (type == NORM_FROBENIUS) { 1610 if (lda>m) { 1611 for (j=0; j<A->cmap->n; j++) { 1612 v = mat->v+j*lda; 1613 for (i=0; i<m; i++) { 1614 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1615 } 1616 } 1617 } else { 1618 #if defined(PETSC_USE_REAL___FP16) 1619 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1620 *nrm = BLASnrm2_(&cnt,v,&one); 1621 } 1622 #else 1623 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1624 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1625 } 1626 } 1627 *nrm = PetscSqrtReal(sum); 1628 #endif 1629 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1630 } else if (type == NORM_1) { 1631 *nrm = 0.0; 1632 for (j=0; j<A->cmap->n; j++) { 1633 v = mat->v + j*mat->lda; 1634 sum = 0.0; 1635 for (i=0; i<A->rmap->n; i++) { 1636 sum += PetscAbsScalar(*v); v++; 1637 } 1638 if (sum > *nrm) *nrm = sum; 1639 } 1640 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1641 } else if (type == NORM_INFINITY) { 1642 *nrm = 0.0; 1643 for (j=0; j<A->rmap->n; j++) { 1644 v = mat->v + j; 1645 sum = 0.0; 1646 for (i=0; i<A->cmap->n; i++) { 1647 sum += PetscAbsScalar(*v); v += mat->lda; 1648 } 1649 if (sum > *nrm) *nrm = sum; 1650 } 1651 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1652 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1653 PetscFunctionReturn(0); 1654 } 1655 1656 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1657 { 1658 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1659 PetscErrorCode ierr; 1660 1661 PetscFunctionBegin; 1662 switch (op) { 1663 case MAT_ROW_ORIENTED: 1664 aij->roworiented = flg; 1665 break; 1666 case MAT_NEW_NONZERO_LOCATIONS: 1667 case MAT_NEW_NONZERO_LOCATION_ERR: 1668 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1669 case MAT_NEW_DIAGONALS: 1670 case MAT_KEEP_NONZERO_PATTERN: 1671 case MAT_IGNORE_OFF_PROC_ENTRIES: 1672 case MAT_USE_HASH_TABLE: 1673 case MAT_IGNORE_ZERO_ENTRIES: 1674 case MAT_IGNORE_LOWER_TRIANGULAR: 1675 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1676 break; 1677 case MAT_SPD: 1678 case MAT_SYMMETRIC: 1679 case MAT_STRUCTURALLY_SYMMETRIC: 1680 case MAT_HERMITIAN: 1681 case MAT_SYMMETRY_ETERNAL: 1682 /* These options are handled directly by MatSetOption() */ 1683 break; 1684 default: 1685 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1686 } 1687 PetscFunctionReturn(0); 1688 } 1689 1690 static PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1691 { 1692 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1693 PetscErrorCode ierr; 1694 PetscInt lda=l->lda,m=A->rmap->n,j; 1695 1696 PetscFunctionBegin; 1697 if (lda>m) { 1698 for (j=0; j<A->cmap->n; j++) { 1699 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1700 } 1701 } else { 1702 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1703 } 1704 PetscFunctionReturn(0); 1705 } 1706 1707 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1708 { 1709 PetscErrorCode ierr; 1710 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1711 PetscInt m = l->lda, n = A->cmap->n, i,j; 1712 PetscScalar *slot,*bb; 1713 const PetscScalar *xx; 1714 1715 PetscFunctionBegin; 1716 #if defined(PETSC_USE_DEBUG) 1717 for (i=0; i<N; i++) { 1718 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1719 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); 1720 } 1721 #endif 1722 1723 /* fix right hand side if needed */ 1724 if (x && b) { 1725 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1726 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1727 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1728 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1729 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1730 } 1731 1732 for (i=0; i<N; i++) { 1733 slot = l->v + rows[i]; 1734 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1735 } 1736 if (diag != 0.0) { 1737 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1738 for (i=0; i<N; i++) { 1739 slot = l->v + (m+1)*rows[i]; 1740 *slot = diag; 1741 } 1742 } 1743 PetscFunctionReturn(0); 1744 } 1745 1746 static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1747 { 1748 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1749 1750 PetscFunctionBegin; 1751 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"); 1752 *array = mat->v; 1753 PetscFunctionReturn(0); 1754 } 1755 1756 static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1757 { 1758 PetscFunctionBegin; 1759 PetscFunctionReturn(0); 1760 } 1761 1762 /*@C 1763 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1764 1765 Logically Collective on Mat 1766 1767 Input Parameter: 1768 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1769 1770 Output Parameter: 1771 . array - pointer to the data 1772 1773 Level: intermediate 1774 1775 .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1776 @*/ 1777 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1778 { 1779 PetscErrorCode ierr; 1780 1781 PetscFunctionBegin; 1782 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1783 PetscFunctionReturn(0); 1784 } 1785 1786 /*@C 1787 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1788 1789 Logically Collective on Mat 1790 1791 Input Parameters: 1792 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1793 . array - pointer to the data 1794 1795 Level: intermediate 1796 1797 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1798 @*/ 1799 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1800 { 1801 PetscErrorCode ierr; 1802 1803 PetscFunctionBegin; 1804 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1805 if (array) *array = NULL; 1806 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1807 PetscFunctionReturn(0); 1808 } 1809 1810 /*@C 1811 MatDenseGetArrayRead - gives access to the array where the data for a SeqDense matrix is stored 1812 1813 Not Collective 1814 1815 Input Parameter: 1816 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1817 1818 Output Parameter: 1819 . array - pointer to the data 1820 1821 Level: intermediate 1822 1823 .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead() 1824 @*/ 1825 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 1826 { 1827 PetscErrorCode ierr; 1828 1829 PetscFunctionBegin; 1830 ierr = PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1831 PetscFunctionReturn(0); 1832 } 1833 1834 /*@C 1835 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1836 1837 Not Collective 1838 1839 Input Parameters: 1840 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1841 . array - pointer to the data 1842 1843 Level: intermediate 1844 1845 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray() 1846 @*/ 1847 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 1848 { 1849 PetscErrorCode ierr; 1850 1851 PetscFunctionBegin; 1852 ierr = PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1853 if (array) *array = NULL; 1854 PetscFunctionReturn(0); 1855 } 1856 1857 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1858 { 1859 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1860 PetscErrorCode ierr; 1861 PetscInt i,j,nrows,ncols; 1862 const PetscInt *irow,*icol; 1863 PetscScalar *av,*bv,*v = mat->v; 1864 Mat newmat; 1865 1866 PetscFunctionBegin; 1867 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1868 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1869 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1870 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1871 1872 /* Check submatrixcall */ 1873 if (scall == MAT_REUSE_MATRIX) { 1874 PetscInt n_cols,n_rows; 1875 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1876 if (n_rows != nrows || n_cols != ncols) { 1877 /* resize the result matrix to match number of requested rows/columns */ 1878 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1879 } 1880 newmat = *B; 1881 } else { 1882 /* Create and fill new matrix */ 1883 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1884 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1885 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1886 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1887 } 1888 1889 /* Now extract the data pointers and do the copy,column at a time */ 1890 bv = ((Mat_SeqDense*)newmat->data)->v; 1891 1892 for (i=0; i<ncols; i++) { 1893 av = v + mat->lda*icol[i]; 1894 for (j=0; j<nrows; j++) *bv++ = av[irow[j]]; 1895 } 1896 1897 /* Assemble the matrices so that the correct flags are set */ 1898 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1899 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1900 1901 /* Free work space */ 1902 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1903 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1904 *B = newmat; 1905 PetscFunctionReturn(0); 1906 } 1907 1908 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1909 { 1910 PetscErrorCode ierr; 1911 PetscInt i; 1912 1913 PetscFunctionBegin; 1914 if (scall == MAT_INITIAL_MATRIX) { 1915 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 1916 } 1917 1918 for (i=0; i<n; i++) { 1919 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1920 } 1921 PetscFunctionReturn(0); 1922 } 1923 1924 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1925 { 1926 PetscFunctionBegin; 1927 PetscFunctionReturn(0); 1928 } 1929 1930 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1931 { 1932 PetscFunctionBegin; 1933 PetscFunctionReturn(0); 1934 } 1935 1936 static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1937 { 1938 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 1939 PetscErrorCode ierr; 1940 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1941 1942 PetscFunctionBegin; 1943 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1944 if (A->ops->copy != B->ops->copy) { 1945 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1946 PetscFunctionReturn(0); 1947 } 1948 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1949 if (lda1>m || lda2>m) { 1950 for (j=0; j<n; j++) { 1951 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1952 } 1953 } else { 1954 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1955 } 1956 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 1957 PetscFunctionReturn(0); 1958 } 1959 1960 static PetscErrorCode MatSetUp_SeqDense(Mat A) 1961 { 1962 PetscErrorCode ierr; 1963 1964 PetscFunctionBegin; 1965 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1966 PetscFunctionReturn(0); 1967 } 1968 1969 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1970 { 1971 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1972 PetscInt i,nz = A->rmap->n*A->cmap->n; 1973 PetscScalar *aa = a->v; 1974 1975 PetscFunctionBegin; 1976 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1977 PetscFunctionReturn(0); 1978 } 1979 1980 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1981 { 1982 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1983 PetscInt i,nz = A->rmap->n*A->cmap->n; 1984 PetscScalar *aa = a->v; 1985 1986 PetscFunctionBegin; 1987 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1988 PetscFunctionReturn(0); 1989 } 1990 1991 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 1992 { 1993 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1994 PetscInt i,nz = A->rmap->n*A->cmap->n; 1995 PetscScalar *aa = a->v; 1996 1997 PetscFunctionBegin; 1998 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1999 PetscFunctionReturn(0); 2000 } 2001 2002 /* ----------------------------------------------------------------*/ 2003 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2004 { 2005 PetscErrorCode ierr; 2006 2007 PetscFunctionBegin; 2008 if (scall == MAT_INITIAL_MATRIX) { 2009 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2010 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2011 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2012 } 2013 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2014 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2015 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2016 PetscFunctionReturn(0); 2017 } 2018 2019 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2020 { 2021 PetscErrorCode ierr; 2022 PetscInt m=A->rmap->n,n=B->cmap->n; 2023 Mat Cmat; 2024 2025 PetscFunctionBegin; 2026 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); 2027 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2028 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2029 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 2030 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2031 2032 *C = Cmat; 2033 PetscFunctionReturn(0); 2034 } 2035 2036 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2037 { 2038 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2039 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2040 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2041 PetscBLASInt m,n,k; 2042 PetscScalar _DOne=1.0,_DZero=0.0; 2043 PetscErrorCode ierr; 2044 PetscBool flg; 2045 2046 PetscFunctionBegin; 2047 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);CHKERRQ(ierr); 2048 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Second matrix must be dense"); 2049 2050 /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/ 2051 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 2052 if (flg) { 2053 C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense; 2054 ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr); 2055 PetscFunctionReturn(0); 2056 } 2057 2058 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);CHKERRQ(ierr); 2059 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense"); 2060 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2061 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2062 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2063 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2064 PetscFunctionReturn(0); 2065 } 2066 2067 PetscErrorCode MatMatTransposeMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2068 { 2069 PetscErrorCode ierr; 2070 2071 PetscFunctionBegin; 2072 if (scall == MAT_INITIAL_MATRIX) { 2073 ierr = PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2074 ierr = MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2075 ierr = PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2076 } 2077 ierr = PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 2078 ierr = MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2079 ierr = PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 2080 PetscFunctionReturn(0); 2081 } 2082 2083 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2084 { 2085 PetscErrorCode ierr; 2086 PetscInt m=A->rmap->n,n=B->rmap->n; 2087 Mat Cmat; 2088 2089 PetscFunctionBegin; 2090 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); 2091 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2092 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2093 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 2094 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2095 2096 Cmat->assembled = PETSC_TRUE; 2097 2098 *C = Cmat; 2099 PetscFunctionReturn(0); 2100 } 2101 2102 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2103 { 2104 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2105 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2106 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2107 PetscBLASInt m,n,k; 2108 PetscScalar _DOne=1.0,_DZero=0.0; 2109 PetscErrorCode ierr; 2110 2111 PetscFunctionBegin; 2112 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 2113 ierr = PetscBLASIntCast(B->rmap->n,&n);CHKERRQ(ierr); 2114 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2115 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2116 PetscFunctionReturn(0); 2117 } 2118 2119 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2120 { 2121 PetscErrorCode ierr; 2122 2123 PetscFunctionBegin; 2124 if (scall == MAT_INITIAL_MATRIX) { 2125 ierr = PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2126 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2127 ierr = PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2128 } 2129 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2130 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2131 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2132 PetscFunctionReturn(0); 2133 } 2134 2135 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2136 { 2137 PetscErrorCode ierr; 2138 PetscInt m=A->cmap->n,n=B->cmap->n; 2139 Mat Cmat; 2140 2141 PetscFunctionBegin; 2142 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); 2143 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2144 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2145 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 2146 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2147 2148 Cmat->assembled = PETSC_TRUE; 2149 2150 *C = Cmat; 2151 PetscFunctionReturn(0); 2152 } 2153 2154 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2155 { 2156 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2157 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2158 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2159 PetscBLASInt m,n,k; 2160 PetscScalar _DOne=1.0,_DZero=0.0; 2161 PetscErrorCode ierr; 2162 2163 PetscFunctionBegin; 2164 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2165 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2166 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 2167 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2168 PetscFunctionReturn(0); 2169 } 2170 2171 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2172 { 2173 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2174 PetscErrorCode ierr; 2175 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2176 PetscScalar *x; 2177 MatScalar *aa = a->v; 2178 2179 PetscFunctionBegin; 2180 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2181 2182 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2183 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2184 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2185 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2186 for (i=0; i<m; i++) { 2187 x[i] = aa[i]; if (idx) idx[i] = 0; 2188 for (j=1; j<n; j++) { 2189 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 2190 } 2191 } 2192 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2193 PetscFunctionReturn(0); 2194 } 2195 2196 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2197 { 2198 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2199 PetscErrorCode ierr; 2200 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2201 PetscScalar *x; 2202 PetscReal atmp; 2203 MatScalar *aa = a->v; 2204 2205 PetscFunctionBegin; 2206 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2207 2208 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2209 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2210 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2211 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2212 for (i=0; i<m; i++) { 2213 x[i] = PetscAbsScalar(aa[i]); 2214 for (j=1; j<n; j++) { 2215 atmp = PetscAbsScalar(aa[i+m*j]); 2216 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2217 } 2218 } 2219 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2220 PetscFunctionReturn(0); 2221 } 2222 2223 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2224 { 2225 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2226 PetscErrorCode ierr; 2227 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2228 PetscScalar *x; 2229 MatScalar *aa = a->v; 2230 2231 PetscFunctionBegin; 2232 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2233 2234 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2235 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2236 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2237 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2238 for (i=0; i<m; i++) { 2239 x[i] = aa[i]; if (idx) idx[i] = 0; 2240 for (j=1; j<n; j++) { 2241 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 2242 } 2243 } 2244 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2245 PetscFunctionReturn(0); 2246 } 2247 2248 static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2249 { 2250 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2251 PetscErrorCode ierr; 2252 PetscScalar *x; 2253 2254 PetscFunctionBegin; 2255 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2256 2257 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2258 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 2259 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2260 PetscFunctionReturn(0); 2261 } 2262 2263 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2264 { 2265 PetscErrorCode ierr; 2266 PetscInt i,j,m,n; 2267 PetscScalar *a; 2268 2269 PetscFunctionBegin; 2270 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2271 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 2272 ierr = MatDenseGetArray(A,&a);CHKERRQ(ierr); 2273 if (type == NORM_2) { 2274 for (i=0; i<n; i++) { 2275 for (j=0; j<m; j++) { 2276 norms[i] += PetscAbsScalar(a[j]*a[j]); 2277 } 2278 a += m; 2279 } 2280 } else if (type == NORM_1) { 2281 for (i=0; i<n; i++) { 2282 for (j=0; j<m; j++) { 2283 norms[i] += PetscAbsScalar(a[j]); 2284 } 2285 a += m; 2286 } 2287 } else if (type == NORM_INFINITY) { 2288 for (i=0; i<n; i++) { 2289 for (j=0; j<m; j++) { 2290 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2291 } 2292 a += m; 2293 } 2294 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2295 ierr = MatDenseRestoreArray(A,&a);CHKERRQ(ierr); 2296 if (type == NORM_2) { 2297 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2298 } 2299 PetscFunctionReturn(0); 2300 } 2301 2302 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2303 { 2304 PetscErrorCode ierr; 2305 PetscScalar *a; 2306 PetscInt m,n,i; 2307 2308 PetscFunctionBegin; 2309 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2310 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2311 for (i=0; i<m*n; i++) { 2312 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 2313 } 2314 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2315 PetscFunctionReturn(0); 2316 } 2317 2318 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2319 { 2320 PetscFunctionBegin; 2321 *missing = PETSC_FALSE; 2322 PetscFunctionReturn(0); 2323 } 2324 2325 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2326 { 2327 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2328 2329 PetscFunctionBegin; 2330 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2331 *vals = a->v+col*a->lda; 2332 PetscFunctionReturn(0); 2333 } 2334 2335 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2336 { 2337 PetscFunctionBegin; 2338 *vals = 0; /* user cannot accidently use the array later */ 2339 PetscFunctionReturn(0); 2340 } 2341 2342 /* -------------------------------------------------------------------*/ 2343 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2344 MatGetRow_SeqDense, 2345 MatRestoreRow_SeqDense, 2346 MatMult_SeqDense, 2347 /* 4*/ MatMultAdd_SeqDense, 2348 MatMultTranspose_SeqDense, 2349 MatMultTransposeAdd_SeqDense, 2350 0, 2351 0, 2352 0, 2353 /* 10*/ 0, 2354 MatLUFactor_SeqDense, 2355 MatCholeskyFactor_SeqDense, 2356 MatSOR_SeqDense, 2357 MatTranspose_SeqDense, 2358 /* 15*/ MatGetInfo_SeqDense, 2359 MatEqual_SeqDense, 2360 MatGetDiagonal_SeqDense, 2361 MatDiagonalScale_SeqDense, 2362 MatNorm_SeqDense, 2363 /* 20*/ MatAssemblyBegin_SeqDense, 2364 MatAssemblyEnd_SeqDense, 2365 MatSetOption_SeqDense, 2366 MatZeroEntries_SeqDense, 2367 /* 24*/ MatZeroRows_SeqDense, 2368 0, 2369 0, 2370 0, 2371 0, 2372 /* 29*/ MatSetUp_SeqDense, 2373 0, 2374 0, 2375 0, 2376 0, 2377 /* 34*/ MatDuplicate_SeqDense, 2378 0, 2379 0, 2380 0, 2381 0, 2382 /* 39*/ MatAXPY_SeqDense, 2383 MatCreateSubMatrices_SeqDense, 2384 0, 2385 MatGetValues_SeqDense, 2386 MatCopy_SeqDense, 2387 /* 44*/ MatGetRowMax_SeqDense, 2388 MatScale_SeqDense, 2389 MatShift_Basic, 2390 0, 2391 MatZeroRowsColumns_SeqDense, 2392 /* 49*/ MatSetRandom_SeqDense, 2393 0, 2394 0, 2395 0, 2396 0, 2397 /* 54*/ 0, 2398 0, 2399 0, 2400 0, 2401 0, 2402 /* 59*/ 0, 2403 MatDestroy_SeqDense, 2404 MatView_SeqDense, 2405 0, 2406 0, 2407 /* 64*/ 0, 2408 0, 2409 0, 2410 0, 2411 0, 2412 /* 69*/ MatGetRowMaxAbs_SeqDense, 2413 0, 2414 0, 2415 0, 2416 0, 2417 /* 74*/ 0, 2418 0, 2419 0, 2420 0, 2421 0, 2422 /* 79*/ 0, 2423 0, 2424 0, 2425 0, 2426 /* 83*/ MatLoad_SeqDense, 2427 0, 2428 MatIsHermitian_SeqDense, 2429 0, 2430 0, 2431 0, 2432 /* 89*/ MatMatMult_SeqDense_SeqDense, 2433 MatMatMultSymbolic_SeqDense_SeqDense, 2434 MatMatMultNumeric_SeqDense_SeqDense, 2435 MatPtAP_SeqDense_SeqDense, 2436 MatPtAPSymbolic_SeqDense_SeqDense, 2437 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense, 2438 MatMatTransposeMult_SeqDense_SeqDense, 2439 MatMatTransposeMultSymbolic_SeqDense_SeqDense, 2440 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2441 0, 2442 /* 99*/ 0, 2443 0, 2444 0, 2445 MatConjugate_SeqDense, 2446 0, 2447 /*104*/ 0, 2448 MatRealPart_SeqDense, 2449 MatImaginaryPart_SeqDense, 2450 0, 2451 0, 2452 /*109*/ 0, 2453 0, 2454 MatGetRowMin_SeqDense, 2455 MatGetColumnVector_SeqDense, 2456 MatMissingDiagonal_SeqDense, 2457 /*114*/ 0, 2458 0, 2459 0, 2460 0, 2461 0, 2462 /*119*/ 0, 2463 0, 2464 0, 2465 0, 2466 0, 2467 /*124*/ 0, 2468 MatGetColumnNorms_SeqDense, 2469 0, 2470 0, 2471 0, 2472 /*129*/ 0, 2473 MatTransposeMatMult_SeqDense_SeqDense, 2474 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2475 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2476 0, 2477 /*134*/ 0, 2478 0, 2479 0, 2480 0, 2481 0, 2482 /*139*/ 0, 2483 0, 2484 0, 2485 0, 2486 0, 2487 /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqDense 2488 }; 2489 2490 /*@C 2491 MatCreateSeqDense - Creates a sequential dense matrix that 2492 is stored in column major order (the usual Fortran 77 manner). Many 2493 of the matrix operations use the BLAS and LAPACK routines. 2494 2495 Collective on MPI_Comm 2496 2497 Input Parameters: 2498 + comm - MPI communicator, set to PETSC_COMM_SELF 2499 . m - number of rows 2500 . n - number of columns 2501 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2502 to control all matrix memory allocation. 2503 2504 Output Parameter: 2505 . A - the matrix 2506 2507 Notes: 2508 The data input variable is intended primarily for Fortran programmers 2509 who wish to allocate their own matrix memory space. Most users should 2510 set data=NULL. 2511 2512 Level: intermediate 2513 2514 .keywords: dense, matrix, LAPACK, BLAS 2515 2516 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2517 @*/ 2518 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2519 { 2520 PetscErrorCode ierr; 2521 2522 PetscFunctionBegin; 2523 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2524 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2525 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2526 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2527 PetscFunctionReturn(0); 2528 } 2529 2530 /*@C 2531 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2532 2533 Collective on MPI_Comm 2534 2535 Input Parameters: 2536 + B - the matrix 2537 - data - the array (or NULL) 2538 2539 Notes: 2540 The data input variable is intended primarily for Fortran programmers 2541 who wish to allocate their own matrix memory space. Most users should 2542 need not call this routine. 2543 2544 Level: intermediate 2545 2546 .keywords: dense, matrix, LAPACK, BLAS 2547 2548 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2549 2550 @*/ 2551 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2552 { 2553 PetscErrorCode ierr; 2554 2555 PetscFunctionBegin; 2556 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2557 PetscFunctionReturn(0); 2558 } 2559 2560 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2561 { 2562 Mat_SeqDense *b; 2563 PetscErrorCode ierr; 2564 2565 PetscFunctionBegin; 2566 B->preallocated = PETSC_TRUE; 2567 2568 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2569 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2570 2571 b = (Mat_SeqDense*)B->data; 2572 b->Mmax = B->rmap->n; 2573 b->Nmax = B->cmap->n; 2574 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2575 2576 ierr = PetscIntMultError(b->lda,b->Nmax,NULL);CHKERRQ(ierr); 2577 if (!data) { /* petsc-allocated storage */ 2578 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2579 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2580 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2581 2582 b->user_alloc = PETSC_FALSE; 2583 } else { /* user-allocated storage */ 2584 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2585 b->v = data; 2586 b->user_alloc = PETSC_TRUE; 2587 } 2588 B->assembled = PETSC_TRUE; 2589 PetscFunctionReturn(0); 2590 } 2591 2592 #if defined(PETSC_HAVE_ELEMENTAL) 2593 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2594 { 2595 Mat mat_elemental; 2596 PetscErrorCode ierr; 2597 PetscScalar *array,*v_colwise; 2598 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2599 2600 PetscFunctionBegin; 2601 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2602 ierr = MatDenseGetArray(A,&array);CHKERRQ(ierr); 2603 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2604 k = 0; 2605 for (j=0; j<N; j++) { 2606 cols[j] = j; 2607 for (i=0; i<M; i++) { 2608 v_colwise[j*M+i] = array[k++]; 2609 } 2610 } 2611 for (i=0; i<M; i++) { 2612 rows[i] = i; 2613 } 2614 ierr = MatDenseRestoreArray(A,&array);CHKERRQ(ierr); 2615 2616 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2617 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2618 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2619 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2620 2621 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2622 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2623 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2624 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2625 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2626 2627 if (reuse == MAT_INPLACE_MATRIX) { 2628 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2629 } else { 2630 *newmat = mat_elemental; 2631 } 2632 PetscFunctionReturn(0); 2633 } 2634 #endif 2635 2636 /*@C 2637 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2638 2639 Input parameter: 2640 + A - the matrix 2641 - lda - the leading dimension 2642 2643 Notes: 2644 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2645 it asserts that the preallocation has a leading dimension (the LDA parameter 2646 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2647 2648 Level: intermediate 2649 2650 .keywords: dense, matrix, LAPACK, BLAS 2651 2652 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2653 2654 @*/ 2655 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2656 { 2657 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2658 2659 PetscFunctionBegin; 2660 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); 2661 b->lda = lda; 2662 b->changelda = PETSC_FALSE; 2663 b->Mmax = PetscMax(b->Mmax,lda); 2664 PetscFunctionReturn(0); 2665 } 2666 2667 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2668 { 2669 PetscErrorCode ierr; 2670 PetscMPIInt size; 2671 2672 PetscFunctionBegin; 2673 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2674 if (size == 1) { 2675 if (scall == MAT_INITIAL_MATRIX) { 2676 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 2677 } else { 2678 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2679 } 2680 } else { 2681 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2682 } 2683 PetscFunctionReturn(0); 2684 } 2685 2686 /*MC 2687 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2688 2689 Options Database Keys: 2690 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2691 2692 Level: beginner 2693 2694 .seealso: MatCreateSeqDense() 2695 2696 M*/ 2697 2698 PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B) 2699 { 2700 Mat_SeqDense *b; 2701 PetscErrorCode ierr; 2702 PetscMPIInt size; 2703 2704 PetscFunctionBegin; 2705 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2706 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2707 2708 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2709 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2710 B->data = (void*)b; 2711 2712 b->roworiented = PETSC_TRUE; 2713 2714 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2715 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2716 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 2717 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 2718 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2719 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2720 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2721 #if defined(PETSC_HAVE_ELEMENTAL) 2722 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 2723 #endif 2724 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2725 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2726 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2727 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2728 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2729 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2730 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2731 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2732 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2733 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijsell_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2734 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2735 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2736 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijsell_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2737 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2738 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2739 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2740 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2741 2742 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2743 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2744 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2745 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2746 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2747 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2748 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijsell_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2749 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2750 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2751 2752 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2753 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2754 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2755 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr); 2756 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr); 2757 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2758 PetscFunctionReturn(0); 2759 } 2760 2761 /*@C 2762 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. 2763 2764 Not Collective 2765 2766 Input Parameter: 2767 + mat - a MATSEQDENSE or MATMPIDENSE matrix 2768 - col - column index 2769 2770 Output Parameter: 2771 . vals - pointer to the data 2772 2773 Level: intermediate 2774 2775 .seealso: MatDenseRestoreColumn() 2776 @*/ 2777 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 2778 { 2779 PetscErrorCode ierr; 2780 2781 PetscFunctionBegin; 2782 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr); 2783 PetscFunctionReturn(0); 2784 } 2785 2786 /*@C 2787 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 2788 2789 Not Collective 2790 2791 Input Parameter: 2792 . mat - a MATSEQDENSE or MATMPIDENSE matrix 2793 2794 Output Parameter: 2795 . vals - pointer to the data 2796 2797 Level: intermediate 2798 2799 .seealso: MatDenseGetColumn() 2800 @*/ 2801 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 2802 { 2803 PetscErrorCode ierr; 2804 2805 PetscFunctionBegin; 2806 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr); 2807 PetscFunctionReturn(0); 2808 } 2809