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