1 #define PETSCMAT_DLL 2 3 /* 4 Defines the basic matrix operations for sequential dense. 5 */ 6 7 #include "../src/mat/impls/dense/seq/dense.h" 8 #include "petscblaslapack.h" 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "MatAXPY_SeqDense" 12 PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str) 13 { 14 Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data; 15 PetscScalar oalpha = alpha; 16 PetscInt j; 17 PetscBLASInt N,m,ldax,lday,one = 1; 18 PetscErrorCode ierr; 19 20 PetscFunctionBegin; 21 N = PetscBLASIntCast(X->rmap->n*X->cmap->n); 22 m = PetscBLASIntCast(X->rmap->n); 23 ldax = PetscBLASIntCast(x->lda); 24 lday = PetscBLASIntCast(y->lda); 25 if (ldax>m || lday>m) { 26 for (j=0; j<X->cmap->n; j++) { 27 BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one); 28 } 29 } else { 30 BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one); 31 } 32 ierr = PetscLogFlops(PetscMax(2*N-1,0));CHKERRQ(ierr); 33 PetscFunctionReturn(0); 34 } 35 36 #undef __FUNCT__ 37 #define __FUNCT__ "MatGetInfo_SeqDense" 38 PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 39 { 40 PetscInt N = A->rmap->n*A->cmap->n; 41 42 PetscFunctionBegin; 43 info->block_size = 1.0; 44 info->nz_allocated = (double)N; 45 info->nz_used = (double)N; 46 info->nz_unneeded = (double)0; 47 info->assemblies = (double)A->num_ass; 48 info->mallocs = 0; 49 info->memory = ((PetscObject)A)->mem; 50 info->fill_ratio_given = 0; 51 info->fill_ratio_needed = 0; 52 info->factor_mallocs = 0; 53 PetscFunctionReturn(0); 54 } 55 56 #undef __FUNCT__ 57 #define __FUNCT__ "MatScale_SeqDense" 58 PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha) 59 { 60 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 61 PetscScalar oalpha = alpha; 62 PetscErrorCode ierr; 63 PetscBLASInt one = 1,j,nz,lda = PetscBLASIntCast(a->lda); 64 65 PetscFunctionBegin; 66 if (lda>A->rmap->n) { 67 nz = PetscBLASIntCast(A->rmap->n); 68 for (j=0; j<A->cmap->n; j++) { 69 BLASscal_(&nz,&oalpha,a->v+j*lda,&one); 70 } 71 } else { 72 nz = PetscBLASIntCast(A->rmap->n*A->cmap->n); 73 BLASscal_(&nz,&oalpha,a->v,&one); 74 } 75 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 76 PetscFunctionReturn(0); 77 } 78 79 #undef __FUNCT__ 80 #define __FUNCT__ "MatIsHermitian_SeqDense" 81 PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 82 { 83 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 84 PetscInt i,j,m = A->rmap->n,N; 85 PetscScalar *v = a->v; 86 87 PetscFunctionBegin; 88 *fl = PETSC_FALSE; 89 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 90 N = a->lda; 91 92 for (i=0; i<m; i++) { 93 for (j=i+1; j<m; j++) { 94 if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) PetscFunctionReturn(0); 95 } 96 } 97 *fl = PETSC_TRUE; 98 PetscFunctionReturn(0); 99 } 100 101 #undef __FUNCT__ 102 #define __FUNCT__ "MatDuplicateNoCreate_SeqDense" 103 PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues) 104 { 105 Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l; 106 PetscErrorCode ierr; 107 PetscInt lda = (PetscInt)mat->lda,j,m; 108 109 PetscFunctionBegin; 110 ierr = MatSeqDenseSetPreallocation(newi,PETSC_NULL);CHKERRQ(ierr); 111 if (cpvalues == MAT_COPY_VALUES) { 112 l = (Mat_SeqDense*)newi->data; 113 if (lda>A->rmap->n) { 114 m = A->rmap->n; 115 for (j=0; j<A->cmap->n; j++) { 116 ierr = PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 117 } 118 } else { 119 ierr = PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 120 } 121 } 122 newi->assembled = PETSC_TRUE; 123 PetscFunctionReturn(0); 124 } 125 126 #undef __FUNCT__ 127 #define __FUNCT__ "MatDuplicate_SeqDense" 128 PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 129 { 130 PetscErrorCode ierr; 131 132 PetscFunctionBegin; 133 ierr = MatCreate(((PetscObject)A)->comm,newmat);CHKERRQ(ierr); 134 ierr = MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 135 ierr = MatSetType(*newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 136 ierr = MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);CHKERRQ(ierr); 137 PetscFunctionReturn(0); 138 } 139 140 141 extern PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*); 142 143 #undef __FUNCT__ 144 #define __FUNCT__ "MatLUFactorNumeric_SeqDense" 145 PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 146 { 147 MatFactorInfo info; 148 PetscErrorCode ierr; 149 150 PetscFunctionBegin; 151 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 152 ierr = MatLUFactor_SeqDense(fact,0,0,&info);CHKERRQ(ierr); 153 PetscFunctionReturn(0); 154 } 155 156 #undef __FUNCT__ 157 #define __FUNCT__ "MatSolve_SeqDense" 158 PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 159 { 160 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 161 PetscErrorCode ierr; 162 PetscScalar *x,*y; 163 PetscBLASInt one = 1,info,m = PetscBLASIntCast(A->rmap->n); 164 165 PetscFunctionBegin; 166 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 167 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 168 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 169 if (A->factortype == MAT_FACTOR_LU) { 170 #if defined(PETSC_MISSING_LAPACK_GETRS) 171 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 172 #else 173 LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 174 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 175 #endif 176 } else if (A->factortype == MAT_FACTOR_CHOLESKY){ 177 #if defined(PETSC_MISSING_LAPACK_POTRS) 178 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 179 #else 180 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 181 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 182 #endif 183 } 184 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 185 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 186 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 187 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 188 PetscFunctionReturn(0); 189 } 190 191 #undef __FUNCT__ 192 #define __FUNCT__ "MatSolveTranspose_SeqDense" 193 PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy) 194 { 195 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 196 PetscErrorCode ierr; 197 PetscScalar *x,*y; 198 PetscBLASInt one = 1,info,m = PetscBLASIntCast(A->rmap->n); 199 200 PetscFunctionBegin; 201 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 202 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 203 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 204 /* assume if pivots exist then use LU; else Cholesky */ 205 if (mat->pivots) { 206 #if defined(PETSC_MISSING_LAPACK_GETRS) 207 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 208 #else 209 LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 210 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 211 #endif 212 } else { 213 #if defined(PETSC_MISSING_LAPACK_POTRS) 214 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 215 #else 216 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 217 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 218 #endif 219 } 220 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 221 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 222 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 223 PetscFunctionReturn(0); 224 } 225 226 #undef __FUNCT__ 227 #define __FUNCT__ "MatSolveAdd_SeqDense" 228 PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 229 { 230 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 231 PetscErrorCode ierr; 232 PetscScalar *x,*y,sone = 1.0; 233 Vec tmp = 0; 234 PetscBLASInt one = 1,info,m = PetscBLASIntCast(A->rmap->n); 235 236 PetscFunctionBegin; 237 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 238 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 239 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 240 if (yy == zz) { 241 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 242 ierr = PetscLogObjectParent(A,tmp);CHKERRQ(ierr); 243 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 244 } 245 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 246 /* assume if pivots exist then use LU; else Cholesky */ 247 if (mat->pivots) { 248 #if defined(PETSC_MISSING_LAPACK_GETRS) 249 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 250 #else 251 LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 252 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 253 #endif 254 } else { 255 #if defined(PETSC_MISSING_LAPACK_POTRS) 256 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 257 #else 258 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 259 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 260 #endif 261 } 262 if (tmp) {ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); ierr = VecDestroy(tmp);CHKERRQ(ierr);} 263 else {ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr);} 264 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 265 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 266 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n);CHKERRQ(ierr); 267 PetscFunctionReturn(0); 268 } 269 270 #undef __FUNCT__ 271 #define __FUNCT__ "MatSolveTransposeAdd_SeqDense" 272 PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 273 { 274 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 275 PetscErrorCode ierr; 276 PetscScalar *x,*y,sone = 1.0; 277 Vec tmp; 278 PetscBLASInt one = 1,info,m = PetscBLASIntCast(A->rmap->n); 279 280 PetscFunctionBegin; 281 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 282 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 283 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 284 if (yy == zz) { 285 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 286 ierr = PetscLogObjectParent(A,tmp);CHKERRQ(ierr); 287 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 288 } 289 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 290 /* assume if pivots exist then use LU; else Cholesky */ 291 if (mat->pivots) { 292 #if defined(PETSC_MISSING_LAPACK_GETRS) 293 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 294 #else 295 LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 296 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 297 #endif 298 } else { 299 #if defined(PETSC_MISSING_LAPACK_POTRS) 300 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 301 #else 302 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 303 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 304 #endif 305 } 306 if (tmp) { 307 ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); 308 ierr = VecDestroy(tmp);CHKERRQ(ierr); 309 } else { 310 ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr); 311 } 312 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 313 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 314 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n);CHKERRQ(ierr); 315 PetscFunctionReturn(0); 316 } 317 318 /* ---------------------------------------------------------------*/ 319 /* COMMENT: I have chosen to hide row permutation in the pivots, 320 rather than put it in the Mat->row slot.*/ 321 #undef __FUNCT__ 322 #define __FUNCT__ "MatLUFactor_SeqDense" 323 PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo) 324 { 325 #if defined(PETSC_MISSING_LAPACK_GETRF) 326 PetscFunctionBegin; 327 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable."); 328 #else 329 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 330 PetscErrorCode ierr; 331 PetscBLASInt n,m,info; 332 333 PetscFunctionBegin; 334 n = PetscBLASIntCast(A->cmap->n); 335 m = PetscBLASIntCast(A->rmap->n); 336 if (!mat->pivots) { 337 ierr = PetscMalloc((A->rmap->n+1)*sizeof(PetscBLASInt),&mat->pivots);CHKERRQ(ierr); 338 ierr = PetscLogObjectMemory(A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 339 } 340 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 341 LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info); 342 if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization"); 343 if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); 344 A->ops->solve = MatSolve_SeqDense; 345 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 346 A->ops->solveadd = MatSolveAdd_SeqDense; 347 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 348 A->factortype = MAT_FACTOR_LU; 349 350 ierr = PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);CHKERRQ(ierr); 351 #endif 352 PetscFunctionReturn(0); 353 } 354 355 #undef __FUNCT__ 356 #define __FUNCT__ "MatCholeskyFactor_SeqDense" 357 PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo) 358 { 359 #if defined(PETSC_MISSING_LAPACK_POTRF) 360 PetscFunctionBegin; 361 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable."); 362 #else 363 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 364 PetscErrorCode ierr; 365 PetscBLASInt info,n = PetscBLASIntCast(A->cmap->n); 366 367 PetscFunctionBegin; 368 ierr = PetscFree(mat->pivots);CHKERRQ(ierr); 369 mat->pivots = 0; 370 371 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 372 LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info); 373 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1); 374 A->ops->solve = MatSolve_SeqDense; 375 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 376 A->ops->solveadd = MatSolveAdd_SeqDense; 377 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 378 A->factortype = MAT_FACTOR_CHOLESKY; 379 ierr = PetscLogFlops((A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 380 #endif 381 PetscFunctionReturn(0); 382 } 383 384 385 #undef __FUNCT__ 386 #define __FUNCT__ "MatCholeskyFactorNumeric_SeqDense" 387 PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 388 { 389 PetscErrorCode ierr; 390 MatFactorInfo info; 391 392 PetscFunctionBegin; 393 info.fill = 1.0; 394 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 395 ierr = MatCholeskyFactor_SeqDense(fact,0,&info);CHKERRQ(ierr); 396 PetscFunctionReturn(0); 397 } 398 399 #undef __FUNCT__ 400 #define __FUNCT__ "MatCholeskyFactorSymbolic_SeqDense" 401 PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info) 402 { 403 PetscFunctionBegin; 404 fact->assembled = PETSC_TRUE; 405 fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense; 406 PetscFunctionReturn(0); 407 } 408 409 #undef __FUNCT__ 410 #define __FUNCT__ "MatLUFactorSymbolic_SeqDense" 411 PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info) 412 { 413 PetscFunctionBegin; 414 fact->assembled = PETSC_TRUE; 415 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense; 416 PetscFunctionReturn(0); 417 } 418 419 EXTERN_C_BEGIN 420 #undef __FUNCT__ 421 #define __FUNCT__ "MatGetFactor_seqdense_petsc" 422 PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact) 423 { 424 PetscErrorCode ierr; 425 426 PetscFunctionBegin; 427 ierr = MatCreate(((PetscObject)A)->comm,fact);CHKERRQ(ierr); 428 ierr = MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 429 ierr = MatSetType(*fact,((PetscObject)A)->type_name);CHKERRQ(ierr); 430 if (ftype == MAT_FACTOR_LU){ 431 (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense; 432 } else { 433 (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense; 434 } 435 (*fact)->factortype = ftype; 436 PetscFunctionReturn(0); 437 } 438 EXTERN_C_END 439 440 /* ------------------------------------------------------------------*/ 441 #undef __FUNCT__ 442 #define __FUNCT__ "MatSOR_SeqDense" 443 PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx) 444 { 445 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 446 PetscScalar *x,*b,*v = mat->v,zero = 0.0,xt; 447 PetscErrorCode ierr; 448 PetscInt m = A->rmap->n,i; 449 #if !defined(PETSC_USE_COMPLEX) 450 PetscBLASInt o = 1,bm = PetscBLASIntCast(m); 451 #endif 452 453 PetscFunctionBegin; 454 if (flag & SOR_ZERO_INITIAL_GUESS) { 455 /* this is a hack fix, should have another version without the second BLASdot */ 456 ierr = VecSet(xx,zero);CHKERRQ(ierr); 457 } 458 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 459 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 460 its = its*lits; 461 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 462 while (its--) { 463 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 464 for (i=0; i<m; i++) { 465 #if defined(PETSC_USE_COMPLEX) 466 /* cannot use BLAS dot for complex because compiler/linker is 467 not happy about returning a double complex */ 468 PetscInt _i; 469 PetscScalar sum = b[i]; 470 for (_i=0; _i<m; _i++) { 471 sum -= PetscConj(v[i+_i*m])*x[_i]; 472 } 473 xt = sum; 474 #else 475 xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o); 476 #endif 477 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 478 } 479 } 480 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 481 for (i=m-1; i>=0; i--) { 482 #if defined(PETSC_USE_COMPLEX) 483 /* cannot use BLAS dot for complex because compiler/linker is 484 not happy about returning a double complex */ 485 PetscInt _i; 486 PetscScalar sum = b[i]; 487 for (_i=0; _i<m; _i++) { 488 sum -= PetscConj(v[i+_i*m])*x[_i]; 489 } 490 xt = sum; 491 #else 492 xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o); 493 #endif 494 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 495 } 496 } 497 } 498 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 499 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 500 PetscFunctionReturn(0); 501 } 502 503 /* -----------------------------------------------------------------*/ 504 #undef __FUNCT__ 505 #define __FUNCT__ "MatMultTranspose_SeqDense" 506 PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 507 { 508 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 509 PetscScalar *v = mat->v,*x,*y; 510 PetscErrorCode ierr; 511 PetscBLASInt m, n,_One=1; 512 PetscScalar _DOne=1.0,_DZero=0.0; 513 514 PetscFunctionBegin; 515 m = PetscBLASIntCast(A->rmap->n); 516 n = PetscBLASIntCast(A->cmap->n); 517 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 518 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 519 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 520 BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One); 521 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 522 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 523 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 524 PetscFunctionReturn(0); 525 } 526 527 #undef __FUNCT__ 528 #define __FUNCT__ "MatMult_SeqDense" 529 PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 530 { 531 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 532 PetscScalar *v = mat->v,*x,*y,_DOne=1.0,_DZero=0.0; 533 PetscErrorCode ierr; 534 PetscBLASInt m, n, _One=1; 535 536 PetscFunctionBegin; 537 m = PetscBLASIntCast(A->rmap->n); 538 n = PetscBLASIntCast(A->cmap->n); 539 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 540 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 541 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 542 BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One); 543 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 544 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 545 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr); 546 PetscFunctionReturn(0); 547 } 548 549 #undef __FUNCT__ 550 #define __FUNCT__ "MatMultAdd_SeqDense" 551 PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 552 { 553 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 554 PetscScalar *v = mat->v,*x,*y,_DOne=1.0; 555 PetscErrorCode ierr; 556 PetscBLASInt m, n, _One=1; 557 558 PetscFunctionBegin; 559 m = PetscBLASIntCast(A->rmap->n); 560 n = PetscBLASIntCast(A->cmap->n); 561 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 562 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 563 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 564 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 565 BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One); 566 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 567 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 568 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 569 PetscFunctionReturn(0); 570 } 571 572 #undef __FUNCT__ 573 #define __FUNCT__ "MatMultTransposeAdd_SeqDense" 574 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 575 { 576 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 577 PetscScalar *v = mat->v,*x,*y; 578 PetscErrorCode ierr; 579 PetscBLASInt m, n, _One=1; 580 PetscScalar _DOne=1.0; 581 582 PetscFunctionBegin; 583 m = PetscBLASIntCast(A->rmap->n); 584 n = PetscBLASIntCast(A->cmap->n); 585 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 586 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 587 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 588 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 589 BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One); 590 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 591 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 592 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 593 PetscFunctionReturn(0); 594 } 595 596 /* -----------------------------------------------------------------*/ 597 #undef __FUNCT__ 598 #define __FUNCT__ "MatGetRow_SeqDense" 599 PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 600 { 601 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 602 PetscScalar *v; 603 PetscErrorCode ierr; 604 PetscInt i; 605 606 PetscFunctionBegin; 607 *ncols = A->cmap->n; 608 if (cols) { 609 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscInt),cols);CHKERRQ(ierr); 610 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 611 } 612 if (vals) { 613 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscScalar),vals);CHKERRQ(ierr); 614 v = mat->v + row; 615 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 616 } 617 PetscFunctionReturn(0); 618 } 619 620 #undef __FUNCT__ 621 #define __FUNCT__ "MatRestoreRow_SeqDense" 622 PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 623 { 624 PetscErrorCode ierr; 625 PetscFunctionBegin; 626 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 627 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 628 PetscFunctionReturn(0); 629 } 630 /* ----------------------------------------------------------------*/ 631 #undef __FUNCT__ 632 #define __FUNCT__ "MatSetValues_SeqDense" 633 PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 634 { 635 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 636 PetscInt i,j,idx=0; 637 638 PetscFunctionBegin; 639 if (v) PetscValidScalarPointer(v,6); 640 if (!mat->roworiented) { 641 if (addv == INSERT_VALUES) { 642 for (j=0; j<n; j++) { 643 if (indexn[j] < 0) {idx += m; continue;} 644 #if defined(PETSC_USE_DEBUG) 645 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); 646 #endif 647 for (i=0; i<m; i++) { 648 if (indexm[i] < 0) {idx++; continue;} 649 #if defined(PETSC_USE_DEBUG) 650 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); 651 #endif 652 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 653 } 654 } 655 } else { 656 for (j=0; j<n; j++) { 657 if (indexn[j] < 0) {idx += m; continue;} 658 #if defined(PETSC_USE_DEBUG) 659 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); 660 #endif 661 for (i=0; i<m; i++) { 662 if (indexm[i] < 0) {idx++; continue;} 663 #if defined(PETSC_USE_DEBUG) 664 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); 665 #endif 666 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 667 } 668 } 669 } 670 } else { 671 if (addv == INSERT_VALUES) { 672 for (i=0; i<m; i++) { 673 if (indexm[i] < 0) { idx += n; continue;} 674 #if defined(PETSC_USE_DEBUG) 675 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); 676 #endif 677 for (j=0; j<n; j++) { 678 if (indexn[j] < 0) { idx++; continue;} 679 #if defined(PETSC_USE_DEBUG) 680 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); 681 #endif 682 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 683 } 684 } 685 } else { 686 for (i=0; i<m; i++) { 687 if (indexm[i] < 0) { idx += n; continue;} 688 #if defined(PETSC_USE_DEBUG) 689 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); 690 #endif 691 for (j=0; j<n; j++) { 692 if (indexn[j] < 0) { idx++; continue;} 693 #if defined(PETSC_USE_DEBUG) 694 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); 695 #endif 696 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 697 } 698 } 699 } 700 } 701 PetscFunctionReturn(0); 702 } 703 704 #undef __FUNCT__ 705 #define __FUNCT__ "MatGetValues_SeqDense" 706 PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 707 { 708 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 709 PetscInt i,j; 710 711 PetscFunctionBegin; 712 /* row-oriented output */ 713 for (i=0; i<m; i++) { 714 if (indexm[i] < 0) {v += n;continue;} 715 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); 716 for (j=0; j<n; j++) { 717 if (indexn[j] < 0) {v++; continue;} 718 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); 719 *v++ = mat->v[indexn[j]*mat->lda + indexm[i]]; 720 } 721 } 722 PetscFunctionReturn(0); 723 } 724 725 /* -----------------------------------------------------------------*/ 726 727 #undef __FUNCT__ 728 #define __FUNCT__ "MatLoad_SeqDense" 729 PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer) 730 { 731 Mat_SeqDense *a; 732 PetscErrorCode ierr; 733 PetscInt *scols,i,j,nz,header[4]; 734 int fd; 735 PetscMPIInt size; 736 PetscInt *rowlengths = 0,M,N,*cols,grows,gcols; 737 PetscScalar *vals,*svals,*v,*w; 738 MPI_Comm comm = ((PetscObject)viewer)->comm; 739 740 PetscFunctionBegin; 741 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 742 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 743 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 744 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 745 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 746 M = header[1]; N = header[2]; nz = header[3]; 747 748 /* set global size if not set already*/ 749 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 750 ierr = MatSetSizes(newmat,M,N,M,N);CHKERRQ(ierr); 751 } else { 752 /* if sizes and type are already set, check if the vector global sizes are correct */ 753 ierr = MatGetSize(newmat,&grows,&gcols);CHKERRQ(ierr); 754 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); 755 } 756 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 757 758 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 759 a = (Mat_SeqDense*)newmat->data; 760 v = a->v; 761 /* Allocate some temp space to read in the values and then flip them 762 from row major to column major */ 763 ierr = PetscMalloc((M*N > 0 ? M*N : 1)*sizeof(PetscScalar),&w);CHKERRQ(ierr); 764 /* read in nonzero values */ 765 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 766 /* now flip the values and store them in the matrix*/ 767 for (j=0; j<N; j++) { 768 for (i=0; i<M; i++) { 769 *v++ =w[i*N+j]; 770 } 771 } 772 ierr = PetscFree(w);CHKERRQ(ierr); 773 } else { 774 /* read row lengths */ 775 ierr = PetscMalloc((M+1)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 776 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 777 778 a = (Mat_SeqDense*)newmat->data; 779 v = a->v; 780 781 /* read column indices and nonzeros */ 782 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&scols);CHKERRQ(ierr); 783 cols = scols; 784 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 785 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&svals);CHKERRQ(ierr); 786 vals = svals; 787 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 788 789 /* insert into matrix */ 790 for (i=0; i<M; i++) { 791 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 792 svals += rowlengths[i]; scols += rowlengths[i]; 793 } 794 ierr = PetscFree(vals);CHKERRQ(ierr); 795 ierr = PetscFree(cols);CHKERRQ(ierr); 796 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 797 } 798 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 799 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 800 801 PetscFunctionReturn(0); 802 } 803 804 #undef __FUNCT__ 805 #define __FUNCT__ "MatView_SeqDense_ASCII" 806 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 807 { 808 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 809 PetscErrorCode ierr; 810 PetscInt i,j; 811 const char *name; 812 PetscScalar *v; 813 PetscViewerFormat format; 814 #if defined(PETSC_USE_COMPLEX) 815 PetscBool allreal = PETSC_TRUE; 816 #endif 817 818 PetscFunctionBegin; 819 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 820 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 821 PetscFunctionReturn(0); /* do nothing for now */ 822 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 823 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 824 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");CHKERRQ(ierr); 825 for (i=0; i<A->rmap->n; i++) { 826 v = a->v + i; 827 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 828 for (j=0; j<A->cmap->n; j++) { 829 #if defined(PETSC_USE_COMPLEX) 830 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 831 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",j,PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 832 } else if (PetscRealPart(*v)) { 833 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,PetscRealPart(*v));CHKERRQ(ierr); 834 } 835 #else 836 if (*v) { 837 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,*v);CHKERRQ(ierr); 838 } 839 #endif 840 v += a->lda; 841 } 842 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 843 } 844 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 845 } else { 846 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 847 #if defined(PETSC_USE_COMPLEX) 848 /* determine if matrix has all real values */ 849 v = a->v; 850 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 851 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;} 852 } 853 #endif 854 if (format == PETSC_VIEWER_ASCII_MATLAB) { 855 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 856 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 857 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 858 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 859 } else { 860 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");CHKERRQ(ierr); 861 } 862 863 for (i=0; i<A->rmap->n; i++) { 864 v = a->v + i; 865 for (j=0; j<A->cmap->n; j++) { 866 #if defined(PETSC_USE_COMPLEX) 867 if (allreal) { 868 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",PetscRealPart(*v));CHKERRQ(ierr); 869 } else { 870 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16e i ",PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 871 } 872 #else 873 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",*v);CHKERRQ(ierr); 874 #endif 875 v += a->lda; 876 } 877 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 878 } 879 if (format == PETSC_VIEWER_ASCII_MATLAB) { 880 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 881 } 882 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 883 } 884 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 885 PetscFunctionReturn(0); 886 } 887 888 #undef __FUNCT__ 889 #define __FUNCT__ "MatView_SeqDense_Binary" 890 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 891 { 892 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 893 PetscErrorCode ierr; 894 int fd; 895 PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n; 896 PetscScalar *v,*anonz,*vals; 897 PetscViewerFormat format; 898 899 PetscFunctionBegin; 900 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 901 902 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 903 if (format == PETSC_VIEWER_NATIVE) { 904 /* store the matrix as a dense matrix */ 905 ierr = PetscMalloc(4*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 906 col_lens[0] = MAT_FILE_CLASSID; 907 col_lens[1] = m; 908 col_lens[2] = n; 909 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 910 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 911 ierr = PetscFree(col_lens);CHKERRQ(ierr); 912 913 /* write out matrix, by rows */ 914 ierr = PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 915 v = a->v; 916 for (j=0; j<n; j++) { 917 for (i=0; i<m; i++) { 918 vals[j + i*n] = *v++; 919 } 920 } 921 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 922 ierr = PetscFree(vals);CHKERRQ(ierr); 923 } else { 924 ierr = PetscMalloc((4+nz)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 925 col_lens[0] = MAT_FILE_CLASSID; 926 col_lens[1] = m; 927 col_lens[2] = n; 928 col_lens[3] = nz; 929 930 /* store lengths of each row and write (including header) to file */ 931 for (i=0; i<m; i++) col_lens[4+i] = n; 932 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 933 934 /* Possibly should write in smaller increments, not whole matrix at once? */ 935 /* store column indices (zero start index) */ 936 ict = 0; 937 for (i=0; i<m; i++) { 938 for (j=0; j<n; j++) col_lens[ict++] = j; 939 } 940 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 941 ierr = PetscFree(col_lens);CHKERRQ(ierr); 942 943 /* store nonzero values */ 944 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);CHKERRQ(ierr); 945 ict = 0; 946 for (i=0; i<m; i++) { 947 v = a->v + i; 948 for (j=0; j<n; j++) { 949 anonz[ict++] = *v; v += a->lda; 950 } 951 } 952 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 953 ierr = PetscFree(anonz);CHKERRQ(ierr); 954 } 955 PetscFunctionReturn(0); 956 } 957 958 #undef __FUNCT__ 959 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 960 PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 961 { 962 Mat A = (Mat) Aa; 963 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 964 PetscErrorCode ierr; 965 PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j; 966 PetscScalar *v = a->v; 967 PetscViewer viewer; 968 PetscDraw popup; 969 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 970 PetscViewerFormat format; 971 972 PetscFunctionBegin; 973 974 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 975 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 976 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 977 978 /* Loop over matrix elements drawing boxes */ 979 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 980 /* Blue for negative and Red for positive */ 981 color = PETSC_DRAW_BLUE; 982 for(j = 0; j < n; j++) { 983 x_l = j; 984 x_r = x_l + 1.0; 985 for(i = 0; i < m; i++) { 986 y_l = m - i - 1.0; 987 y_r = y_l + 1.0; 988 #if defined(PETSC_USE_COMPLEX) 989 if (PetscRealPart(v[j*m+i]) > 0.) { 990 color = PETSC_DRAW_RED; 991 } else if (PetscRealPart(v[j*m+i]) < 0.) { 992 color = PETSC_DRAW_BLUE; 993 } else { 994 continue; 995 } 996 #else 997 if (v[j*m+i] > 0.) { 998 color = PETSC_DRAW_RED; 999 } else if (v[j*m+i] < 0.) { 1000 color = PETSC_DRAW_BLUE; 1001 } else { 1002 continue; 1003 } 1004 #endif 1005 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1006 } 1007 } 1008 } else { 1009 /* use contour shading to indicate magnitude of values */ 1010 /* first determine max of all nonzero values */ 1011 for(i = 0; i < m*n; i++) { 1012 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1013 } 1014 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 1015 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1016 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 1017 for(j = 0; j < n; j++) { 1018 x_l = j; 1019 x_r = x_l + 1.0; 1020 for(i = 0; i < m; i++) { 1021 y_l = m - i - 1.0; 1022 y_r = y_l + 1.0; 1023 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 1024 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1025 } 1026 } 1027 } 1028 PetscFunctionReturn(0); 1029 } 1030 1031 #undef __FUNCT__ 1032 #define __FUNCT__ "MatView_SeqDense_Draw" 1033 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1034 { 1035 PetscDraw draw; 1036 PetscBool isnull; 1037 PetscReal xr,yr,xl,yl,h,w; 1038 PetscErrorCode ierr; 1039 1040 PetscFunctionBegin; 1041 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1042 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1043 if (isnull) PetscFunctionReturn(0); 1044 1045 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1046 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1047 xr += w; yr += h; xl = -w; yl = -h; 1048 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1049 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1050 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 1051 PetscFunctionReturn(0); 1052 } 1053 1054 #undef __FUNCT__ 1055 #define __FUNCT__ "MatView_SeqDense" 1056 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1057 { 1058 PetscErrorCode ierr; 1059 PetscBool iascii,isbinary,isdraw; 1060 1061 PetscFunctionBegin; 1062 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1063 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1064 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1065 1066 if (iascii) { 1067 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1068 } else if (isbinary) { 1069 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1070 } else if (isdraw) { 1071 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1072 } else { 1073 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name); 1074 } 1075 PetscFunctionReturn(0); 1076 } 1077 1078 #undef __FUNCT__ 1079 #define __FUNCT__ "MatDestroy_SeqDense" 1080 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1081 { 1082 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1083 PetscErrorCode ierr; 1084 1085 PetscFunctionBegin; 1086 #if defined(PETSC_USE_LOG) 1087 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1088 #endif 1089 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1090 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1091 ierr = PetscFree(l);CHKERRQ(ierr); 1092 1093 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1094 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 1095 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1096 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1097 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1098 PetscFunctionReturn(0); 1099 } 1100 1101 #undef __FUNCT__ 1102 #define __FUNCT__ "MatTranspose_SeqDense" 1103 PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1104 { 1105 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1106 PetscErrorCode ierr; 1107 PetscInt k,j,m,n,M; 1108 PetscScalar *v,tmp; 1109 1110 PetscFunctionBegin; 1111 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1112 if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */ 1113 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1114 else { 1115 for (j=0; j<m; j++) { 1116 for (k=0; k<j; k++) { 1117 tmp = v[j + k*M]; 1118 v[j + k*M] = v[k + j*M]; 1119 v[k + j*M] = tmp; 1120 } 1121 } 1122 } 1123 } else { /* out-of-place transpose */ 1124 Mat tmat; 1125 Mat_SeqDense *tmatd; 1126 PetscScalar *v2; 1127 1128 if (reuse == MAT_INITIAL_MATRIX) { 1129 ierr = MatCreate(((PetscObject)A)->comm,&tmat);CHKERRQ(ierr); 1130 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1131 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1132 ierr = MatSeqDenseSetPreallocation(tmat,PETSC_NULL);CHKERRQ(ierr); 1133 } else { 1134 tmat = *matout; 1135 } 1136 tmatd = (Mat_SeqDense*)tmat->data; 1137 v = mat->v; v2 = tmatd->v; 1138 for (j=0; j<n; j++) { 1139 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1140 } 1141 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1142 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1143 *matout = tmat; 1144 } 1145 PetscFunctionReturn(0); 1146 } 1147 1148 #undef __FUNCT__ 1149 #define __FUNCT__ "MatEqual_SeqDense" 1150 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1151 { 1152 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1153 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1154 PetscInt i,j; 1155 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1156 1157 PetscFunctionBegin; 1158 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1159 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1160 for (i=0; i<A1->rmap->n; i++) { 1161 v1 = mat1->v+i; v2 = mat2->v+i; 1162 for (j=0; j<A1->cmap->n; j++) { 1163 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1164 v1 += mat1->lda; v2 += mat2->lda; 1165 } 1166 } 1167 *flg = PETSC_TRUE; 1168 PetscFunctionReturn(0); 1169 } 1170 1171 #undef __FUNCT__ 1172 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1173 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1174 { 1175 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1176 PetscErrorCode ierr; 1177 PetscInt i,n,len; 1178 PetscScalar *x,zero = 0.0; 1179 1180 PetscFunctionBegin; 1181 ierr = VecSet(v,zero);CHKERRQ(ierr); 1182 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1183 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1184 len = PetscMin(A->rmap->n,A->cmap->n); 1185 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1186 for (i=0; i<len; i++) { 1187 x[i] = mat->v[i*mat->lda + i]; 1188 } 1189 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1190 PetscFunctionReturn(0); 1191 } 1192 1193 #undef __FUNCT__ 1194 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1195 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1196 { 1197 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1198 PetscScalar *l,*r,x,*v; 1199 PetscErrorCode ierr; 1200 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1201 1202 PetscFunctionBegin; 1203 if (ll) { 1204 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1205 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1206 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1207 for (i=0; i<m; i++) { 1208 x = l[i]; 1209 v = mat->v + i; 1210 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1211 } 1212 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1213 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1214 } 1215 if (rr) { 1216 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1217 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1218 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1219 for (i=0; i<n; i++) { 1220 x = r[i]; 1221 v = mat->v + i*m; 1222 for (j=0; j<m; j++) { (*v++) *= x;} 1223 } 1224 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1225 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1226 } 1227 PetscFunctionReturn(0); 1228 } 1229 1230 #undef __FUNCT__ 1231 #define __FUNCT__ "MatNorm_SeqDense" 1232 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1233 { 1234 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1235 PetscScalar *v = mat->v; 1236 PetscReal sum = 0.0; 1237 PetscInt lda=mat->lda,m=A->rmap->n,i,j; 1238 PetscErrorCode ierr; 1239 1240 PetscFunctionBegin; 1241 if (type == NORM_FROBENIUS) { 1242 if (lda>m) { 1243 for (j=0; j<A->cmap->n; j++) { 1244 v = mat->v+j*lda; 1245 for (i=0; i<m; i++) { 1246 #if defined(PETSC_USE_COMPLEX) 1247 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1248 #else 1249 sum += (*v)*(*v); v++; 1250 #endif 1251 } 1252 } 1253 } else { 1254 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1255 #if defined(PETSC_USE_COMPLEX) 1256 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1257 #else 1258 sum += (*v)*(*v); v++; 1259 #endif 1260 } 1261 } 1262 *nrm = sqrt(sum); 1263 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1264 } else if (type == NORM_1) { 1265 *nrm = 0.0; 1266 for (j=0; j<A->cmap->n; j++) { 1267 v = mat->v + j*mat->lda; 1268 sum = 0.0; 1269 for (i=0; i<A->rmap->n; i++) { 1270 sum += PetscAbsScalar(*v); v++; 1271 } 1272 if (sum > *nrm) *nrm = sum; 1273 } 1274 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1275 } else if (type == NORM_INFINITY) { 1276 *nrm = 0.0; 1277 for (j=0; j<A->rmap->n; j++) { 1278 v = mat->v + j; 1279 sum = 0.0; 1280 for (i=0; i<A->cmap->n; i++) { 1281 sum += PetscAbsScalar(*v); v += mat->lda; 1282 } 1283 if (sum > *nrm) *nrm = sum; 1284 } 1285 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1286 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1287 PetscFunctionReturn(0); 1288 } 1289 1290 #undef __FUNCT__ 1291 #define __FUNCT__ "MatSetOption_SeqDense" 1292 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1293 { 1294 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1295 PetscErrorCode ierr; 1296 1297 PetscFunctionBegin; 1298 switch (op) { 1299 case MAT_ROW_ORIENTED: 1300 aij->roworiented = flg; 1301 break; 1302 case MAT_NEW_NONZERO_LOCATIONS: 1303 case MAT_NEW_NONZERO_LOCATION_ERR: 1304 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1305 case MAT_NEW_DIAGONALS: 1306 case MAT_IGNORE_OFF_PROC_ENTRIES: 1307 case MAT_USE_HASH_TABLE: 1308 case MAT_SYMMETRIC: 1309 case MAT_STRUCTURALLY_SYMMETRIC: 1310 case MAT_HERMITIAN: 1311 case MAT_SYMMETRY_ETERNAL: 1312 case MAT_IGNORE_LOWER_TRIANGULAR: 1313 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1314 break; 1315 default: 1316 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1317 } 1318 PetscFunctionReturn(0); 1319 } 1320 1321 #undef __FUNCT__ 1322 #define __FUNCT__ "MatZeroEntries_SeqDense" 1323 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1324 { 1325 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1326 PetscErrorCode ierr; 1327 PetscInt lda=l->lda,m=A->rmap->n,j; 1328 1329 PetscFunctionBegin; 1330 if (lda>m) { 1331 for (j=0; j<A->cmap->n; j++) { 1332 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1333 } 1334 } else { 1335 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1336 } 1337 PetscFunctionReturn(0); 1338 } 1339 1340 #undef __FUNCT__ 1341 #define __FUNCT__ "MatZeroRows_SeqDense" 1342 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1343 { 1344 PetscErrorCode ierr; 1345 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1346 PetscInt n = A->cmap->n,i,j; 1347 PetscScalar *slot,*bb; 1348 const PetscScalar *xx; 1349 1350 PetscFunctionBegin; 1351 /* fix right hand side if needed */ 1352 if (x && b) { 1353 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1354 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1355 for (i=0; i<N; i++) { 1356 bb[rows[i]] = diag*xx[rows[i]]; 1357 } 1358 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1359 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1360 } 1361 1362 for (i=0; i<N; i++) { 1363 slot = l->v + rows[i]; 1364 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1365 } 1366 if (diag != 0.0) { 1367 for (i=0; i<N; i++) { 1368 slot = l->v + (n+1)*rows[i]; 1369 *slot = diag; 1370 } 1371 } 1372 PetscFunctionReturn(0); 1373 } 1374 1375 #undef __FUNCT__ 1376 #define __FUNCT__ "MatGetArray_SeqDense" 1377 PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[]) 1378 { 1379 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1380 1381 PetscFunctionBegin; 1382 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"); 1383 *array = mat->v; 1384 PetscFunctionReturn(0); 1385 } 1386 1387 #undef __FUNCT__ 1388 #define __FUNCT__ "MatRestoreArray_SeqDense" 1389 PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1390 { 1391 PetscFunctionBegin; 1392 *array = 0; /* user cannot accidently use the array later */ 1393 PetscFunctionReturn(0); 1394 } 1395 1396 #undef __FUNCT__ 1397 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1398 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1399 { 1400 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1401 PetscErrorCode ierr; 1402 PetscInt i,j,nrows,ncols; 1403 const PetscInt *irow,*icol; 1404 PetscScalar *av,*bv,*v = mat->v; 1405 Mat newmat; 1406 1407 PetscFunctionBegin; 1408 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1409 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1410 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1411 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1412 1413 /* Check submatrixcall */ 1414 if (scall == MAT_REUSE_MATRIX) { 1415 PetscInt n_cols,n_rows; 1416 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1417 if (n_rows != nrows || n_cols != ncols) { 1418 /* resize the result result matrix to match number of requested rows/columns */ 1419 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1420 } 1421 newmat = *B; 1422 } else { 1423 /* Create and fill new matrix */ 1424 ierr = MatCreate(((PetscObject)A)->comm,&newmat);CHKERRQ(ierr); 1425 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1426 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1427 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1428 } 1429 1430 /* Now extract the data pointers and do the copy,column at a time */ 1431 bv = ((Mat_SeqDense*)newmat->data)->v; 1432 1433 for (i=0; i<ncols; i++) { 1434 av = v + mat->lda*icol[i]; 1435 for (j=0; j<nrows; j++) { 1436 *bv++ = av[irow[j]]; 1437 } 1438 } 1439 1440 /* Assemble the matrices so that the correct flags are set */ 1441 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1442 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1443 1444 /* Free work space */ 1445 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1446 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1447 *B = newmat; 1448 PetscFunctionReturn(0); 1449 } 1450 1451 #undef __FUNCT__ 1452 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1453 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1454 { 1455 PetscErrorCode ierr; 1456 PetscInt i; 1457 1458 PetscFunctionBegin; 1459 if (scall == MAT_INITIAL_MATRIX) { 1460 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1461 } 1462 1463 for (i=0; i<n; i++) { 1464 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1465 } 1466 PetscFunctionReturn(0); 1467 } 1468 1469 #undef __FUNCT__ 1470 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1471 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1472 { 1473 PetscFunctionBegin; 1474 PetscFunctionReturn(0); 1475 } 1476 1477 #undef __FUNCT__ 1478 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1479 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1480 { 1481 PetscFunctionBegin; 1482 PetscFunctionReturn(0); 1483 } 1484 1485 #undef __FUNCT__ 1486 #define __FUNCT__ "MatCopy_SeqDense" 1487 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1488 { 1489 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1490 PetscErrorCode ierr; 1491 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1492 1493 PetscFunctionBegin; 1494 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1495 if (A->ops->copy != B->ops->copy) { 1496 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1497 PetscFunctionReturn(0); 1498 } 1499 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1500 if (lda1>m || lda2>m) { 1501 for (j=0; j<n; j++) { 1502 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1503 } 1504 } else { 1505 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1506 } 1507 PetscFunctionReturn(0); 1508 } 1509 1510 #undef __FUNCT__ 1511 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1512 PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A) 1513 { 1514 PetscErrorCode ierr; 1515 1516 PetscFunctionBegin; 1517 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1518 PetscFunctionReturn(0); 1519 } 1520 1521 #undef __FUNCT__ 1522 #define __FUNCT__ "MatSetSizes_SeqDense" 1523 PetscErrorCode MatSetSizes_SeqDense(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) 1524 { 1525 PetscFunctionBegin; 1526 /* this will not be called before lda, Mmax, and Nmax have been set */ 1527 m = PetscMax(m,M); 1528 n = PetscMax(n,N); 1529 1530 /* if (m > a->Mmax) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot yet resize number rows of dense matrix larger then its initial size %d, requested %d",a->lda,(int)m); 1531 if (n > a->Nmax) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot yet resize number columns of dense matrix larger then its initial size %d, requested %d",a->Nmax,(int)n); 1532 */ 1533 A->rmap->n = A->rmap->N = m; 1534 A->cmap->n = A->cmap->N = n; 1535 PetscFunctionReturn(0); 1536 } 1537 1538 #undef __FUNCT__ 1539 #define __FUNCT__ "MatConjugate_SeqDense" 1540 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1541 { 1542 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1543 PetscInt i,nz = A->rmap->n*A->cmap->n; 1544 PetscScalar *aa = a->v; 1545 1546 PetscFunctionBegin; 1547 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1548 PetscFunctionReturn(0); 1549 } 1550 1551 #undef __FUNCT__ 1552 #define __FUNCT__ "MatRealPart_SeqDense" 1553 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1554 { 1555 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1556 PetscInt i,nz = A->rmap->n*A->cmap->n; 1557 PetscScalar *aa = a->v; 1558 1559 PetscFunctionBegin; 1560 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1561 PetscFunctionReturn(0); 1562 } 1563 1564 #undef __FUNCT__ 1565 #define __FUNCT__ "MatImaginaryPart_SeqDense" 1566 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 1567 { 1568 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1569 PetscInt i,nz = A->rmap->n*A->cmap->n; 1570 PetscScalar *aa = a->v; 1571 1572 PetscFunctionBegin; 1573 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1574 PetscFunctionReturn(0); 1575 } 1576 1577 /* ----------------------------------------------------------------*/ 1578 #undef __FUNCT__ 1579 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1580 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1581 { 1582 PetscErrorCode ierr; 1583 1584 PetscFunctionBegin; 1585 if (scall == MAT_INITIAL_MATRIX){ 1586 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1587 } 1588 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1589 PetscFunctionReturn(0); 1590 } 1591 1592 #undef __FUNCT__ 1593 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1594 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1595 { 1596 PetscErrorCode ierr; 1597 PetscInt m=A->rmap->n,n=B->cmap->n; 1598 Mat Cmat; 1599 1600 PetscFunctionBegin; 1601 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); 1602 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1603 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1604 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1605 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1606 Cmat->assembled = PETSC_TRUE; 1607 *C = Cmat; 1608 PetscFunctionReturn(0); 1609 } 1610 1611 #undef __FUNCT__ 1612 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1613 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1614 { 1615 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1616 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1617 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1618 PetscBLASInt m,n,k; 1619 PetscScalar _DOne=1.0,_DZero=0.0; 1620 1621 PetscFunctionBegin; 1622 m = PetscBLASIntCast(A->rmap->n); 1623 n = PetscBLASIntCast(B->cmap->n); 1624 k = PetscBLASIntCast(A->cmap->n); 1625 BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1626 PetscFunctionReturn(0); 1627 } 1628 1629 #undef __FUNCT__ 1630 #define __FUNCT__ "MatMatMultTranspose_SeqDense_SeqDense" 1631 PetscErrorCode MatMatMultTranspose_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1632 { 1633 PetscErrorCode ierr; 1634 1635 PetscFunctionBegin; 1636 if (scall == MAT_INITIAL_MATRIX){ 1637 ierr = MatMatMultTransposeSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1638 } 1639 ierr = MatMatMultTransposeNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1640 PetscFunctionReturn(0); 1641 } 1642 1643 #undef __FUNCT__ 1644 #define __FUNCT__ "MatMatMultTransposeSymbolic_SeqDense_SeqDense" 1645 PetscErrorCode MatMatMultTransposeSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1646 { 1647 PetscErrorCode ierr; 1648 PetscInt m=A->cmap->n,n=B->cmap->n; 1649 Mat Cmat; 1650 1651 PetscFunctionBegin; 1652 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); 1653 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1654 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1655 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1656 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1657 Cmat->assembled = PETSC_TRUE; 1658 *C = Cmat; 1659 PetscFunctionReturn(0); 1660 } 1661 1662 #undef __FUNCT__ 1663 #define __FUNCT__ "MatMatMultTransposeNumeric_SeqDense_SeqDense" 1664 PetscErrorCode MatMatMultTransposeNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1665 { 1666 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1667 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1668 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1669 PetscBLASInt m,n,k; 1670 PetscScalar _DOne=1.0,_DZero=0.0; 1671 1672 PetscFunctionBegin; 1673 m = PetscBLASIntCast(A->cmap->n); 1674 n = PetscBLASIntCast(B->cmap->n); 1675 k = PetscBLASIntCast(A->rmap->n); 1676 /* 1677 Note the m and n arguments below are the number rows and columns of A', not A! 1678 */ 1679 BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1680 PetscFunctionReturn(0); 1681 } 1682 1683 #undef __FUNCT__ 1684 #define __FUNCT__ "MatGetRowMax_SeqDense" 1685 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1686 { 1687 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1688 PetscErrorCode ierr; 1689 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1690 PetscScalar *x; 1691 MatScalar *aa = a->v; 1692 1693 PetscFunctionBegin; 1694 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1695 1696 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1697 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1698 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1699 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1700 for (i=0; i<m; i++) { 1701 x[i] = aa[i]; if (idx) idx[i] = 0; 1702 for (j=1; j<n; j++){ 1703 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1704 } 1705 } 1706 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1707 PetscFunctionReturn(0); 1708 } 1709 1710 #undef __FUNCT__ 1711 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1712 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1713 { 1714 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1715 PetscErrorCode ierr; 1716 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1717 PetscScalar *x; 1718 PetscReal atmp; 1719 MatScalar *aa = a->v; 1720 1721 PetscFunctionBegin; 1722 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1723 1724 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1725 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1726 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1727 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1728 for (i=0; i<m; i++) { 1729 x[i] = PetscAbsScalar(aa[i]); 1730 for (j=1; j<n; j++){ 1731 atmp = PetscAbsScalar(aa[i+m*j]); 1732 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1733 } 1734 } 1735 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1736 PetscFunctionReturn(0); 1737 } 1738 1739 #undef __FUNCT__ 1740 #define __FUNCT__ "MatGetRowMin_SeqDense" 1741 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1742 { 1743 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1744 PetscErrorCode ierr; 1745 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1746 PetscScalar *x; 1747 MatScalar *aa = a->v; 1748 1749 PetscFunctionBegin; 1750 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1751 1752 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1753 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1754 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1755 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1756 for (i=0; i<m; i++) { 1757 x[i] = aa[i]; if (idx) idx[i] = 0; 1758 for (j=1; j<n; j++){ 1759 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1760 } 1761 } 1762 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1763 PetscFunctionReturn(0); 1764 } 1765 1766 #undef __FUNCT__ 1767 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1768 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1769 { 1770 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1771 PetscErrorCode ierr; 1772 PetscScalar *x; 1773 1774 PetscFunctionBegin; 1775 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1776 1777 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1778 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1779 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1780 PetscFunctionReturn(0); 1781 } 1782 1783 /* -------------------------------------------------------------------*/ 1784 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1785 MatGetRow_SeqDense, 1786 MatRestoreRow_SeqDense, 1787 MatMult_SeqDense, 1788 /* 4*/ MatMultAdd_SeqDense, 1789 MatMultTranspose_SeqDense, 1790 MatMultTransposeAdd_SeqDense, 1791 0, 1792 0, 1793 0, 1794 /*10*/ 0, 1795 MatLUFactor_SeqDense, 1796 MatCholeskyFactor_SeqDense, 1797 MatSOR_SeqDense, 1798 MatTranspose_SeqDense, 1799 /*15*/ MatGetInfo_SeqDense, 1800 MatEqual_SeqDense, 1801 MatGetDiagonal_SeqDense, 1802 MatDiagonalScale_SeqDense, 1803 MatNorm_SeqDense, 1804 /*20*/ MatAssemblyBegin_SeqDense, 1805 MatAssemblyEnd_SeqDense, 1806 MatSetOption_SeqDense, 1807 MatZeroEntries_SeqDense, 1808 /*24*/ MatZeroRows_SeqDense, 1809 0, 1810 0, 1811 0, 1812 0, 1813 /*29*/ MatSetUpPreallocation_SeqDense, 1814 0, 1815 0, 1816 MatGetArray_SeqDense, 1817 MatRestoreArray_SeqDense, 1818 /*34*/ MatDuplicate_SeqDense, 1819 0, 1820 0, 1821 0, 1822 0, 1823 /*39*/ MatAXPY_SeqDense, 1824 MatGetSubMatrices_SeqDense, 1825 0, 1826 MatGetValues_SeqDense, 1827 MatCopy_SeqDense, 1828 /*44*/ MatGetRowMax_SeqDense, 1829 MatScale_SeqDense, 1830 0, 1831 0, 1832 0, 1833 /*49*/ 0, 1834 0, 1835 0, 1836 0, 1837 0, 1838 /*54*/ 0, 1839 0, 1840 0, 1841 0, 1842 0, 1843 /*59*/ 0, 1844 MatDestroy_SeqDense, 1845 MatView_SeqDense, 1846 0, 1847 0, 1848 /*64*/ 0, 1849 0, 1850 0, 1851 0, 1852 0, 1853 /*69*/ MatGetRowMaxAbs_SeqDense, 1854 0, 1855 0, 1856 0, 1857 0, 1858 /*74*/ 0, 1859 0, 1860 0, 1861 0, 1862 0, 1863 /*79*/ 0, 1864 0, 1865 0, 1866 0, 1867 /*83*/ MatLoad_SeqDense, 1868 0, 1869 MatIsHermitian_SeqDense, 1870 0, 1871 0, 1872 0, 1873 /*89*/ MatMatMult_SeqDense_SeqDense, 1874 MatMatMultSymbolic_SeqDense_SeqDense, 1875 MatMatMultNumeric_SeqDense_SeqDense, 1876 0, 1877 0, 1878 /*94*/ 0, 1879 MatMatMultTranspose_SeqDense_SeqDense, 1880 MatMatMultTransposeSymbolic_SeqDense_SeqDense, 1881 MatMatMultTransposeNumeric_SeqDense_SeqDense, 1882 0, 1883 /*99*/ 0, 1884 0, 1885 0, 1886 MatConjugate_SeqDense, 1887 MatSetSizes_SeqDense, 1888 /*104*/0, 1889 MatRealPart_SeqDense, 1890 MatImaginaryPart_SeqDense, 1891 0, 1892 0, 1893 /*109*/0, 1894 0, 1895 MatGetRowMin_SeqDense, 1896 MatGetColumnVector_SeqDense, 1897 0, 1898 /*114*/0, 1899 0, 1900 0, 1901 0, 1902 0, 1903 /*119*/0, 1904 0, 1905 0, 1906 0 1907 }; 1908 1909 #undef __FUNCT__ 1910 #define __FUNCT__ "MatCreateSeqDense" 1911 /*@C 1912 MatCreateSeqDense - Creates a sequential dense matrix that 1913 is stored in column major order (the usual Fortran 77 manner). Many 1914 of the matrix operations use the BLAS and LAPACK routines. 1915 1916 Collective on MPI_Comm 1917 1918 Input Parameters: 1919 + comm - MPI communicator, set to PETSC_COMM_SELF 1920 . m - number of rows 1921 . n - number of columns 1922 - data - optional location of matrix data in column major order. Set data=PETSC_NULL for PETSc 1923 to control all matrix memory allocation. 1924 1925 Output Parameter: 1926 . A - the matrix 1927 1928 Notes: 1929 The data input variable is intended primarily for Fortran programmers 1930 who wish to allocate their own matrix memory space. Most users should 1931 set data=PETSC_NULL. 1932 1933 Level: intermediate 1934 1935 .keywords: dense, matrix, LAPACK, BLAS 1936 1937 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1938 @*/ 1939 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 1940 { 1941 PetscErrorCode ierr; 1942 1943 PetscFunctionBegin; 1944 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1945 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1946 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1947 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1948 PetscFunctionReturn(0); 1949 } 1950 1951 #undef __FUNCT__ 1952 #define __FUNCT__ "MatSeqDenseSetPreallocation" 1953 /*@C 1954 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1955 1956 Collective on MPI_Comm 1957 1958 Input Parameters: 1959 + A - the matrix 1960 - data - the array (or PETSC_NULL) 1961 1962 Notes: 1963 The data input variable is intended primarily for Fortran programmers 1964 who wish to allocate their own matrix memory space. Most users should 1965 need not call this routine. 1966 1967 Level: intermediate 1968 1969 .keywords: dense, matrix, LAPACK, BLAS 1970 1971 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues(), MatSeqDenseSetLDA() 1972 1973 @*/ 1974 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 1975 { 1976 PetscErrorCode ierr; 1977 1978 PetscFunctionBegin; 1979 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 1980 PetscFunctionReturn(0); 1981 } 1982 1983 EXTERN_C_BEGIN 1984 #undef __FUNCT__ 1985 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 1986 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1987 { 1988 Mat_SeqDense *b; 1989 PetscErrorCode ierr; 1990 1991 PetscFunctionBegin; 1992 B->preallocated = PETSC_TRUE; 1993 1994 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 1995 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 1996 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1997 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1998 1999 b = (Mat_SeqDense*)B->data; 2000 b->Mmax = B->rmap->n; 2001 b->Nmax = B->cmap->n; 2002 if(b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2003 2004 if (!data) { /* petsc-allocated storage */ 2005 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2006 ierr = PetscMalloc(b->lda*b->Nmax*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 2007 ierr = PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2008 ierr = PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2009 b->user_alloc = PETSC_FALSE; 2010 } else { /* user-allocated storage */ 2011 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2012 b->v = data; 2013 b->user_alloc = PETSC_TRUE; 2014 } 2015 B->assembled = PETSC_TRUE; 2016 PetscFunctionReturn(0); 2017 } 2018 EXTERN_C_END 2019 2020 #undef __FUNCT__ 2021 #define __FUNCT__ "MatSeqDenseSetLDA" 2022 /*@C 2023 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2024 2025 Input parameter: 2026 + A - the matrix 2027 - lda - the leading dimension 2028 2029 Notes: 2030 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2031 it asserts that the preallocation has a leading dimension (the LDA parameter 2032 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2033 2034 Level: intermediate 2035 2036 .keywords: dense, matrix, LAPACK, BLAS 2037 2038 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2039 2040 @*/ 2041 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetLDA(Mat B,PetscInt lda) 2042 { 2043 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2044 2045 PetscFunctionBegin; 2046 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); 2047 b->lda = lda; 2048 b->changelda = PETSC_FALSE; 2049 b->Mmax = PetscMax(b->Mmax,lda); 2050 PetscFunctionReturn(0); 2051 } 2052 2053 /*MC 2054 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2055 2056 Options Database Keys: 2057 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2058 2059 Level: beginner 2060 2061 .seealso: MatCreateSeqDense() 2062 2063 M*/ 2064 2065 EXTERN_C_BEGIN 2066 #undef __FUNCT__ 2067 #define __FUNCT__ "MatCreate_SeqDense" 2068 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqDense(Mat B) 2069 { 2070 Mat_SeqDense *b; 2071 PetscErrorCode ierr; 2072 PetscMPIInt size; 2073 2074 PetscFunctionBegin; 2075 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 2076 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2077 2078 ierr = PetscNewLog(B,Mat_SeqDense,&b);CHKERRQ(ierr); 2079 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2080 B->mapping = 0; 2081 B->data = (void*)b; 2082 2083 b->pivots = 0; 2084 b->roworiented = PETSC_TRUE; 2085 b->v = 0; 2086 b->changelda = PETSC_FALSE; 2087 2088 2089 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 2090 "MatGetFactor_seqdense_petsc", 2091 MatGetFactor_seqdense_petsc);CHKERRQ(ierr); 2092 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 2093 "MatSeqDenseSetPreallocation_SeqDense", 2094 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2095 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 2096 "MatMatMult_SeqAIJ_SeqDense", 2097 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2098 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C", 2099 "MatMatMultSymbolic_SeqAIJ_SeqDense", 2100 MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2101 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C", 2102 "MatMatMultNumeric_SeqAIJ_SeqDense", 2103 MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2104 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2105 PetscFunctionReturn(0); 2106 } 2107 EXTERN_C_END 2108