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,PetscTruth *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 PetscTruth allreal = PETSC_TRUE; 816 #endif 817 818 PetscFunctionBegin; 819 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 820 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 821 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 822 PetscFunctionReturn(0); /* do nothing for now */ 823 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 824 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);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 } 860 861 for (i=0; i<A->rmap->n; i++) { 862 v = a->v + i; 863 for (j=0; j<A->cmap->n; j++) { 864 #if defined(PETSC_USE_COMPLEX) 865 if (allreal) { 866 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",PetscRealPart(*v));CHKERRQ(ierr); 867 } else { 868 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16e i ",PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 869 } 870 #else 871 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",*v);CHKERRQ(ierr); 872 #endif 873 v += a->lda; 874 } 875 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 876 } 877 if (format == PETSC_VIEWER_ASCII_MATLAB) { 878 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 879 } 880 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 881 } 882 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 883 PetscFunctionReturn(0); 884 } 885 886 #undef __FUNCT__ 887 #define __FUNCT__ "MatView_SeqDense_Binary" 888 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 889 { 890 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 891 PetscErrorCode ierr; 892 int fd; 893 PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n; 894 PetscScalar *v,*anonz; 895 896 PetscFunctionBegin; 897 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 898 899 ierr = PetscMalloc((4+nz)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 900 col_lens[0] = MAT_FILE_CLASSID; 901 col_lens[1] = m; 902 col_lens[2] = n; 903 col_lens[3] = nz; 904 905 /* store lengths of each row and write (including header) to file */ 906 for (i=0; i<m; i++) col_lens[4+i] = n; 907 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 908 909 /* Possibly should write in smaller increments, not whole matrix at once? */ 910 /* store column indices (zero start index) */ 911 ict = 0; 912 for (i=0; i<m; i++) { 913 for (j=0; j<n; j++) col_lens[ict++] = j; 914 } 915 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 916 ierr = PetscFree(col_lens);CHKERRQ(ierr); 917 918 /* store nonzero values */ 919 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);CHKERRQ(ierr); 920 ict = 0; 921 for (i=0; i<m; i++) { 922 v = a->v + i; 923 for (j=0; j<n; j++) { 924 anonz[ict++] = *v; v += a->lda; 925 } 926 } 927 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 928 ierr = PetscFree(anonz);CHKERRQ(ierr); 929 930 PetscFunctionReturn(0); 931 } 932 933 #undef __FUNCT__ 934 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 935 PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 936 { 937 Mat A = (Mat) Aa; 938 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 939 PetscErrorCode ierr; 940 PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j; 941 PetscScalar *v = a->v; 942 PetscViewer viewer; 943 PetscDraw popup; 944 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 945 PetscViewerFormat format; 946 947 PetscFunctionBegin; 948 949 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 950 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 951 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 952 953 /* Loop over matrix elements drawing boxes */ 954 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 955 /* Blue for negative and Red for positive */ 956 color = PETSC_DRAW_BLUE; 957 for(j = 0; j < n; j++) { 958 x_l = j; 959 x_r = x_l + 1.0; 960 for(i = 0; i < m; i++) { 961 y_l = m - i - 1.0; 962 y_r = y_l + 1.0; 963 #if defined(PETSC_USE_COMPLEX) 964 if (PetscRealPart(v[j*m+i]) > 0.) { 965 color = PETSC_DRAW_RED; 966 } else if (PetscRealPart(v[j*m+i]) < 0.) { 967 color = PETSC_DRAW_BLUE; 968 } else { 969 continue; 970 } 971 #else 972 if (v[j*m+i] > 0.) { 973 color = PETSC_DRAW_RED; 974 } else if (v[j*m+i] < 0.) { 975 color = PETSC_DRAW_BLUE; 976 } else { 977 continue; 978 } 979 #endif 980 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 981 } 982 } 983 } else { 984 /* use contour shading to indicate magnitude of values */ 985 /* first determine max of all nonzero values */ 986 for(i = 0; i < m*n; i++) { 987 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 988 } 989 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 990 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 991 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 992 for(j = 0; j < n; j++) { 993 x_l = j; 994 x_r = x_l + 1.0; 995 for(i = 0; i < m; i++) { 996 y_l = m - i - 1.0; 997 y_r = y_l + 1.0; 998 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 999 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1000 } 1001 } 1002 } 1003 PetscFunctionReturn(0); 1004 } 1005 1006 #undef __FUNCT__ 1007 #define __FUNCT__ "MatView_SeqDense_Draw" 1008 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1009 { 1010 PetscDraw draw; 1011 PetscTruth isnull; 1012 PetscReal xr,yr,xl,yl,h,w; 1013 PetscErrorCode ierr; 1014 1015 PetscFunctionBegin; 1016 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1017 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1018 if (isnull) PetscFunctionReturn(0); 1019 1020 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1021 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1022 xr += w; yr += h; xl = -w; yl = -h; 1023 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1024 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1025 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 1026 PetscFunctionReturn(0); 1027 } 1028 1029 #undef __FUNCT__ 1030 #define __FUNCT__ "MatView_SeqDense" 1031 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1032 { 1033 PetscErrorCode ierr; 1034 PetscTruth iascii,isbinary,isdraw; 1035 1036 PetscFunctionBegin; 1037 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1038 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1039 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1040 1041 if (iascii) { 1042 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1043 } else if (isbinary) { 1044 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1045 } else if (isdraw) { 1046 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1047 } else { 1048 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name); 1049 } 1050 PetscFunctionReturn(0); 1051 } 1052 1053 #undef __FUNCT__ 1054 #define __FUNCT__ "MatDestroy_SeqDense" 1055 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1056 { 1057 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1058 PetscErrorCode ierr; 1059 1060 PetscFunctionBegin; 1061 #if defined(PETSC_USE_LOG) 1062 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1063 #endif 1064 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1065 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1066 ierr = PetscFree(l);CHKERRQ(ierr); 1067 1068 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1069 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 1070 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1071 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1072 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1073 PetscFunctionReturn(0); 1074 } 1075 1076 #undef __FUNCT__ 1077 #define __FUNCT__ "MatTranspose_SeqDense" 1078 PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1079 { 1080 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1081 PetscErrorCode ierr; 1082 PetscInt k,j,m,n,M; 1083 PetscScalar *v,tmp; 1084 1085 PetscFunctionBegin; 1086 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1087 if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */ 1088 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1089 else { 1090 for (j=0; j<m; j++) { 1091 for (k=0; k<j; k++) { 1092 tmp = v[j + k*M]; 1093 v[j + k*M] = v[k + j*M]; 1094 v[k + j*M] = tmp; 1095 } 1096 } 1097 } 1098 } else { /* out-of-place transpose */ 1099 Mat tmat; 1100 Mat_SeqDense *tmatd; 1101 PetscScalar *v2; 1102 1103 if (reuse == MAT_INITIAL_MATRIX) { 1104 ierr = MatCreate(((PetscObject)A)->comm,&tmat);CHKERRQ(ierr); 1105 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1106 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1107 ierr = MatSeqDenseSetPreallocation(tmat,PETSC_NULL);CHKERRQ(ierr); 1108 } else { 1109 tmat = *matout; 1110 } 1111 tmatd = (Mat_SeqDense*)tmat->data; 1112 v = mat->v; v2 = tmatd->v; 1113 for (j=0; j<n; j++) { 1114 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1115 } 1116 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1117 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1118 *matout = tmat; 1119 } 1120 PetscFunctionReturn(0); 1121 } 1122 1123 #undef __FUNCT__ 1124 #define __FUNCT__ "MatEqual_SeqDense" 1125 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg) 1126 { 1127 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1128 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1129 PetscInt i,j; 1130 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1131 1132 PetscFunctionBegin; 1133 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1134 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1135 for (i=0; i<A1->rmap->n; i++) { 1136 v1 = mat1->v+i; v2 = mat2->v+i; 1137 for (j=0; j<A1->cmap->n; j++) { 1138 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1139 v1 += mat1->lda; v2 += mat2->lda; 1140 } 1141 } 1142 *flg = PETSC_TRUE; 1143 PetscFunctionReturn(0); 1144 } 1145 1146 #undef __FUNCT__ 1147 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1148 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1149 { 1150 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1151 PetscErrorCode ierr; 1152 PetscInt i,n,len; 1153 PetscScalar *x,zero = 0.0; 1154 1155 PetscFunctionBegin; 1156 ierr = VecSet(v,zero);CHKERRQ(ierr); 1157 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1158 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1159 len = PetscMin(A->rmap->n,A->cmap->n); 1160 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1161 for (i=0; i<len; i++) { 1162 x[i] = mat->v[i*mat->lda + i]; 1163 } 1164 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1165 PetscFunctionReturn(0); 1166 } 1167 1168 #undef __FUNCT__ 1169 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1170 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1171 { 1172 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1173 PetscScalar *l,*r,x,*v; 1174 PetscErrorCode ierr; 1175 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1176 1177 PetscFunctionBegin; 1178 if (ll) { 1179 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1180 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1181 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1182 for (i=0; i<m; i++) { 1183 x = l[i]; 1184 v = mat->v + i; 1185 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1186 } 1187 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1188 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1189 } 1190 if (rr) { 1191 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1192 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1193 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1194 for (i=0; i<n; i++) { 1195 x = r[i]; 1196 v = mat->v + i*m; 1197 for (j=0; j<m; j++) { (*v++) *= x;} 1198 } 1199 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1200 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1201 } 1202 PetscFunctionReturn(0); 1203 } 1204 1205 #undef __FUNCT__ 1206 #define __FUNCT__ "MatNorm_SeqDense" 1207 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1208 { 1209 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1210 PetscScalar *v = mat->v; 1211 PetscReal sum = 0.0; 1212 PetscInt lda=mat->lda,m=A->rmap->n,i,j; 1213 PetscErrorCode ierr; 1214 1215 PetscFunctionBegin; 1216 if (type == NORM_FROBENIUS) { 1217 if (lda>m) { 1218 for (j=0; j<A->cmap->n; j++) { 1219 v = mat->v+j*lda; 1220 for (i=0; i<m; i++) { 1221 #if defined(PETSC_USE_COMPLEX) 1222 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1223 #else 1224 sum += (*v)*(*v); v++; 1225 #endif 1226 } 1227 } 1228 } else { 1229 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1230 #if defined(PETSC_USE_COMPLEX) 1231 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1232 #else 1233 sum += (*v)*(*v); v++; 1234 #endif 1235 } 1236 } 1237 *nrm = sqrt(sum); 1238 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1239 } else if (type == NORM_1) { 1240 *nrm = 0.0; 1241 for (j=0; j<A->cmap->n; j++) { 1242 v = mat->v + j*mat->lda; 1243 sum = 0.0; 1244 for (i=0; i<A->rmap->n; i++) { 1245 sum += PetscAbsScalar(*v); v++; 1246 } 1247 if (sum > *nrm) *nrm = sum; 1248 } 1249 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1250 } else if (type == NORM_INFINITY) { 1251 *nrm = 0.0; 1252 for (j=0; j<A->rmap->n; j++) { 1253 v = mat->v + j; 1254 sum = 0.0; 1255 for (i=0; i<A->cmap->n; i++) { 1256 sum += PetscAbsScalar(*v); v += mat->lda; 1257 } 1258 if (sum > *nrm) *nrm = sum; 1259 } 1260 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1261 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1262 PetscFunctionReturn(0); 1263 } 1264 1265 #undef __FUNCT__ 1266 #define __FUNCT__ "MatSetOption_SeqDense" 1267 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscTruth flg) 1268 { 1269 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1270 PetscErrorCode ierr; 1271 1272 PetscFunctionBegin; 1273 switch (op) { 1274 case MAT_ROW_ORIENTED: 1275 aij->roworiented = flg; 1276 break; 1277 case MAT_NEW_NONZERO_LOCATIONS: 1278 case MAT_NEW_NONZERO_LOCATION_ERR: 1279 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1280 case MAT_NEW_DIAGONALS: 1281 case MAT_IGNORE_OFF_PROC_ENTRIES: 1282 case MAT_USE_HASH_TABLE: 1283 case MAT_SYMMETRIC: 1284 case MAT_STRUCTURALLY_SYMMETRIC: 1285 case MAT_HERMITIAN: 1286 case MAT_SYMMETRY_ETERNAL: 1287 case MAT_IGNORE_LOWER_TRIANGULAR: 1288 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1289 break; 1290 default: 1291 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1292 } 1293 PetscFunctionReturn(0); 1294 } 1295 1296 #undef __FUNCT__ 1297 #define __FUNCT__ "MatZeroEntries_SeqDense" 1298 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1299 { 1300 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1301 PetscErrorCode ierr; 1302 PetscInt lda=l->lda,m=A->rmap->n,j; 1303 1304 PetscFunctionBegin; 1305 if (lda>m) { 1306 for (j=0; j<A->cmap->n; j++) { 1307 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1308 } 1309 } else { 1310 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1311 } 1312 PetscFunctionReturn(0); 1313 } 1314 1315 #undef __FUNCT__ 1316 #define __FUNCT__ "MatZeroRows_SeqDense" 1317 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1318 { 1319 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1320 PetscInt n = A->cmap->n,i,j; 1321 PetscScalar *slot; 1322 1323 PetscFunctionBegin; 1324 for (i=0; i<N; i++) { 1325 slot = l->v + rows[i]; 1326 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1327 } 1328 if (diag != 0.0) { 1329 for (i=0; i<N; i++) { 1330 slot = l->v + (n+1)*rows[i]; 1331 *slot = diag; 1332 } 1333 } 1334 PetscFunctionReturn(0); 1335 } 1336 1337 #undef __FUNCT__ 1338 #define __FUNCT__ "MatGetArray_SeqDense" 1339 PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[]) 1340 { 1341 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1342 1343 PetscFunctionBegin; 1344 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"); 1345 *array = mat->v; 1346 PetscFunctionReturn(0); 1347 } 1348 1349 #undef __FUNCT__ 1350 #define __FUNCT__ "MatRestoreArray_SeqDense" 1351 PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1352 { 1353 PetscFunctionBegin; 1354 *array = 0; /* user cannot accidently use the array later */ 1355 PetscFunctionReturn(0); 1356 } 1357 1358 #undef __FUNCT__ 1359 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1360 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1361 { 1362 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1363 PetscErrorCode ierr; 1364 PetscInt i,j,nrows,ncols; 1365 const PetscInt *irow,*icol; 1366 PetscScalar *av,*bv,*v = mat->v; 1367 Mat newmat; 1368 1369 PetscFunctionBegin; 1370 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1371 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1372 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1373 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1374 1375 /* Check submatrixcall */ 1376 if (scall == MAT_REUSE_MATRIX) { 1377 PetscInt n_cols,n_rows; 1378 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1379 if (n_rows != nrows || n_cols != ncols) { 1380 /* resize the result result matrix to match number of requested rows/columns */ 1381 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1382 } 1383 newmat = *B; 1384 } else { 1385 /* Create and fill new matrix */ 1386 ierr = MatCreate(((PetscObject)A)->comm,&newmat);CHKERRQ(ierr); 1387 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1388 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1389 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1390 } 1391 1392 /* Now extract the data pointers and do the copy,column at a time */ 1393 bv = ((Mat_SeqDense*)newmat->data)->v; 1394 1395 for (i=0; i<ncols; i++) { 1396 av = v + mat->lda*icol[i]; 1397 for (j=0; j<nrows; j++) { 1398 *bv++ = av[irow[j]]; 1399 } 1400 } 1401 1402 /* Assemble the matrices so that the correct flags are set */ 1403 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1404 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1405 1406 /* Free work space */ 1407 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1408 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1409 *B = newmat; 1410 PetscFunctionReturn(0); 1411 } 1412 1413 #undef __FUNCT__ 1414 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1415 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1416 { 1417 PetscErrorCode ierr; 1418 PetscInt i; 1419 1420 PetscFunctionBegin; 1421 if (scall == MAT_INITIAL_MATRIX) { 1422 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1423 } 1424 1425 for (i=0; i<n; i++) { 1426 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1427 } 1428 PetscFunctionReturn(0); 1429 } 1430 1431 #undef __FUNCT__ 1432 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1433 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1434 { 1435 PetscFunctionBegin; 1436 PetscFunctionReturn(0); 1437 } 1438 1439 #undef __FUNCT__ 1440 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1441 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1442 { 1443 PetscFunctionBegin; 1444 PetscFunctionReturn(0); 1445 } 1446 1447 #undef __FUNCT__ 1448 #define __FUNCT__ "MatCopy_SeqDense" 1449 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1450 { 1451 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1452 PetscErrorCode ierr; 1453 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1454 1455 PetscFunctionBegin; 1456 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1457 if (A->ops->copy != B->ops->copy) { 1458 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1459 PetscFunctionReturn(0); 1460 } 1461 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1462 if (lda1>m || lda2>m) { 1463 for (j=0; j<n; j++) { 1464 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1465 } 1466 } else { 1467 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1468 } 1469 PetscFunctionReturn(0); 1470 } 1471 1472 #undef __FUNCT__ 1473 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1474 PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A) 1475 { 1476 PetscErrorCode ierr; 1477 1478 PetscFunctionBegin; 1479 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1480 PetscFunctionReturn(0); 1481 } 1482 1483 #undef __FUNCT__ 1484 #define __FUNCT__ "MatSetSizes_SeqDense" 1485 PetscErrorCode MatSetSizes_SeqDense(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) 1486 { 1487 PetscFunctionBegin; 1488 /* this will not be called before lda, Mmax, and Nmax have been set */ 1489 m = PetscMax(m,M); 1490 n = PetscMax(n,N); 1491 1492 /* 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); 1493 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); 1494 */ 1495 A->rmap->n = A->rmap->N = m; 1496 A->cmap->n = A->cmap->N = n; 1497 PetscFunctionReturn(0); 1498 } 1499 1500 #undef __FUNCT__ 1501 #define __FUNCT__ "MatConjugate_SeqDense" 1502 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1503 { 1504 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1505 PetscInt i,nz = A->rmap->n*A->cmap->n; 1506 PetscScalar *aa = a->v; 1507 1508 PetscFunctionBegin; 1509 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1510 PetscFunctionReturn(0); 1511 } 1512 1513 #undef __FUNCT__ 1514 #define __FUNCT__ "MatRealPart_SeqDense" 1515 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1516 { 1517 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1518 PetscInt i,nz = A->rmap->n*A->cmap->n; 1519 PetscScalar *aa = a->v; 1520 1521 PetscFunctionBegin; 1522 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1523 PetscFunctionReturn(0); 1524 } 1525 1526 #undef __FUNCT__ 1527 #define __FUNCT__ "MatImaginaryPart_SeqDense" 1528 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 1529 { 1530 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1531 PetscInt i,nz = A->rmap->n*A->cmap->n; 1532 PetscScalar *aa = a->v; 1533 1534 PetscFunctionBegin; 1535 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1536 PetscFunctionReturn(0); 1537 } 1538 1539 /* ----------------------------------------------------------------*/ 1540 #undef __FUNCT__ 1541 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1542 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1543 { 1544 PetscErrorCode ierr; 1545 1546 PetscFunctionBegin; 1547 if (scall == MAT_INITIAL_MATRIX){ 1548 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1549 } 1550 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1551 PetscFunctionReturn(0); 1552 } 1553 1554 #undef __FUNCT__ 1555 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1556 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1557 { 1558 PetscErrorCode ierr; 1559 PetscInt m=A->rmap->n,n=B->cmap->n; 1560 Mat Cmat; 1561 1562 PetscFunctionBegin; 1563 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); 1564 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1565 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1566 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1567 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1568 Cmat->assembled = PETSC_TRUE; 1569 *C = Cmat; 1570 PetscFunctionReturn(0); 1571 } 1572 1573 #undef __FUNCT__ 1574 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1575 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1576 { 1577 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1578 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1579 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1580 PetscBLASInt m,n,k; 1581 PetscScalar _DOne=1.0,_DZero=0.0; 1582 1583 PetscFunctionBegin; 1584 m = PetscBLASIntCast(A->rmap->n); 1585 n = PetscBLASIntCast(B->cmap->n); 1586 k = PetscBLASIntCast(A->cmap->n); 1587 BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1588 PetscFunctionReturn(0); 1589 } 1590 1591 #undef __FUNCT__ 1592 #define __FUNCT__ "MatMatMultTranspose_SeqDense_SeqDense" 1593 PetscErrorCode MatMatMultTranspose_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1594 { 1595 PetscErrorCode ierr; 1596 1597 PetscFunctionBegin; 1598 if (scall == MAT_INITIAL_MATRIX){ 1599 ierr = MatMatMultTransposeSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1600 } 1601 ierr = MatMatMultTransposeNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1602 PetscFunctionReturn(0); 1603 } 1604 1605 #undef __FUNCT__ 1606 #define __FUNCT__ "MatMatMultTransposeSymbolic_SeqDense_SeqDense" 1607 PetscErrorCode MatMatMultTransposeSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1608 { 1609 PetscErrorCode ierr; 1610 PetscInt m=A->cmap->n,n=B->cmap->n; 1611 Mat Cmat; 1612 1613 PetscFunctionBegin; 1614 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); 1615 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1616 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1617 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1618 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1619 Cmat->assembled = PETSC_TRUE; 1620 *C = Cmat; 1621 PetscFunctionReturn(0); 1622 } 1623 1624 #undef __FUNCT__ 1625 #define __FUNCT__ "MatMatMultTransposeNumeric_SeqDense_SeqDense" 1626 PetscErrorCode MatMatMultTransposeNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1627 { 1628 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1629 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1630 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1631 PetscBLASInt m,n,k; 1632 PetscScalar _DOne=1.0,_DZero=0.0; 1633 1634 PetscFunctionBegin; 1635 m = PetscBLASIntCast(A->cmap->n); 1636 n = PetscBLASIntCast(B->cmap->n); 1637 k = PetscBLASIntCast(A->rmap->n); 1638 /* 1639 Note the m and n arguments below are the number rows and columns of A', not A! 1640 */ 1641 BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1642 PetscFunctionReturn(0); 1643 } 1644 1645 #undef __FUNCT__ 1646 #define __FUNCT__ "MatGetRowMax_SeqDense" 1647 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1648 { 1649 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1650 PetscErrorCode ierr; 1651 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1652 PetscScalar *x; 1653 MatScalar *aa = a->v; 1654 1655 PetscFunctionBegin; 1656 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1657 1658 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1659 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1660 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1661 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1662 for (i=0; i<m; i++) { 1663 x[i] = aa[i]; if (idx) idx[i] = 0; 1664 for (j=1; j<n; j++){ 1665 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1666 } 1667 } 1668 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1669 PetscFunctionReturn(0); 1670 } 1671 1672 #undef __FUNCT__ 1673 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1674 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1675 { 1676 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1677 PetscErrorCode ierr; 1678 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1679 PetscScalar *x; 1680 PetscReal atmp; 1681 MatScalar *aa = a->v; 1682 1683 PetscFunctionBegin; 1684 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1685 1686 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1687 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1688 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1689 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1690 for (i=0; i<m; i++) { 1691 x[i] = PetscAbsScalar(aa[i]); 1692 for (j=1; j<n; j++){ 1693 atmp = PetscAbsScalar(aa[i+m*j]); 1694 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1695 } 1696 } 1697 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1698 PetscFunctionReturn(0); 1699 } 1700 1701 #undef __FUNCT__ 1702 #define __FUNCT__ "MatGetRowMin_SeqDense" 1703 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1704 { 1705 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1706 PetscErrorCode ierr; 1707 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1708 PetscScalar *x; 1709 MatScalar *aa = a->v; 1710 1711 PetscFunctionBegin; 1712 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1713 1714 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1715 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1716 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1717 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1718 for (i=0; i<m; i++) { 1719 x[i] = aa[i]; if (idx) idx[i] = 0; 1720 for (j=1; j<n; j++){ 1721 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1722 } 1723 } 1724 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1725 PetscFunctionReturn(0); 1726 } 1727 1728 #undef __FUNCT__ 1729 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1730 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1731 { 1732 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1733 PetscErrorCode ierr; 1734 PetscScalar *x; 1735 1736 PetscFunctionBegin; 1737 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1738 1739 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1740 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1741 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1742 PetscFunctionReturn(0); 1743 } 1744 1745 /* -------------------------------------------------------------------*/ 1746 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1747 MatGetRow_SeqDense, 1748 MatRestoreRow_SeqDense, 1749 MatMult_SeqDense, 1750 /* 4*/ MatMultAdd_SeqDense, 1751 MatMultTranspose_SeqDense, 1752 MatMultTransposeAdd_SeqDense, 1753 0, 1754 0, 1755 0, 1756 /*10*/ 0, 1757 MatLUFactor_SeqDense, 1758 MatCholeskyFactor_SeqDense, 1759 MatSOR_SeqDense, 1760 MatTranspose_SeqDense, 1761 /*15*/ MatGetInfo_SeqDense, 1762 MatEqual_SeqDense, 1763 MatGetDiagonal_SeqDense, 1764 MatDiagonalScale_SeqDense, 1765 MatNorm_SeqDense, 1766 /*20*/ MatAssemblyBegin_SeqDense, 1767 MatAssemblyEnd_SeqDense, 1768 MatSetOption_SeqDense, 1769 MatZeroEntries_SeqDense, 1770 /*24*/ MatZeroRows_SeqDense, 1771 0, 1772 0, 1773 0, 1774 0, 1775 /*29*/ MatSetUpPreallocation_SeqDense, 1776 0, 1777 0, 1778 MatGetArray_SeqDense, 1779 MatRestoreArray_SeqDense, 1780 /*34*/ MatDuplicate_SeqDense, 1781 0, 1782 0, 1783 0, 1784 0, 1785 /*39*/ MatAXPY_SeqDense, 1786 MatGetSubMatrices_SeqDense, 1787 0, 1788 MatGetValues_SeqDense, 1789 MatCopy_SeqDense, 1790 /*44*/ MatGetRowMax_SeqDense, 1791 MatScale_SeqDense, 1792 0, 1793 0, 1794 0, 1795 /*49*/ 0, 1796 0, 1797 0, 1798 0, 1799 0, 1800 /*54*/ 0, 1801 0, 1802 0, 1803 0, 1804 0, 1805 /*59*/ 0, 1806 MatDestroy_SeqDense, 1807 MatView_SeqDense, 1808 0, 1809 0, 1810 /*64*/ 0, 1811 0, 1812 0, 1813 0, 1814 0, 1815 /*69*/ MatGetRowMaxAbs_SeqDense, 1816 0, 1817 0, 1818 0, 1819 0, 1820 /*74*/ 0, 1821 0, 1822 0, 1823 0, 1824 0, 1825 /*79*/ 0, 1826 0, 1827 0, 1828 0, 1829 /*83*/ MatLoad_SeqDense, 1830 0, 1831 MatIsHermitian_SeqDense, 1832 0, 1833 0, 1834 0, 1835 /*89*/ MatMatMult_SeqDense_SeqDense, 1836 MatMatMultSymbolic_SeqDense_SeqDense, 1837 MatMatMultNumeric_SeqDense_SeqDense, 1838 0, 1839 0, 1840 /*94*/ 0, 1841 MatMatMultTranspose_SeqDense_SeqDense, 1842 MatMatMultTransposeSymbolic_SeqDense_SeqDense, 1843 MatMatMultTransposeNumeric_SeqDense_SeqDense, 1844 0, 1845 /*99*/ 0, 1846 0, 1847 0, 1848 MatConjugate_SeqDense, 1849 MatSetSizes_SeqDense, 1850 /*104*/0, 1851 MatRealPart_SeqDense, 1852 MatImaginaryPart_SeqDense, 1853 0, 1854 0, 1855 /*109*/0, 1856 0, 1857 MatGetRowMin_SeqDense, 1858 MatGetColumnVector_SeqDense, 1859 0, 1860 /*114*/0, 1861 0, 1862 0, 1863 0, 1864 0, 1865 /*119*/0, 1866 0, 1867 0, 1868 0 1869 }; 1870 1871 #undef __FUNCT__ 1872 #define __FUNCT__ "MatCreateSeqDense" 1873 /*@C 1874 MatCreateSeqDense - Creates a sequential dense matrix that 1875 is stored in column major order (the usual Fortran 77 manner). Many 1876 of the matrix operations use the BLAS and LAPACK routines. 1877 1878 Collective on MPI_Comm 1879 1880 Input Parameters: 1881 + comm - MPI communicator, set to PETSC_COMM_SELF 1882 . m - number of rows 1883 . n - number of columns 1884 - data - optional location of matrix data in column major order. Set data=PETSC_NULL for PETSc 1885 to control all matrix memory allocation. 1886 1887 Output Parameter: 1888 . A - the matrix 1889 1890 Notes: 1891 The data input variable is intended primarily for Fortran programmers 1892 who wish to allocate their own matrix memory space. Most users should 1893 set data=PETSC_NULL. 1894 1895 Level: intermediate 1896 1897 .keywords: dense, matrix, LAPACK, BLAS 1898 1899 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1900 @*/ 1901 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 1902 { 1903 PetscErrorCode ierr; 1904 1905 PetscFunctionBegin; 1906 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1907 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1908 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1909 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1910 PetscFunctionReturn(0); 1911 } 1912 1913 #undef __FUNCT__ 1914 #define __FUNCT__ "MatSeqDenseSetPreallocation" 1915 /*@C 1916 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1917 1918 Collective on MPI_Comm 1919 1920 Input Parameters: 1921 + A - the matrix 1922 - data - the array (or PETSC_NULL) 1923 1924 Notes: 1925 The data input variable is intended primarily for Fortran programmers 1926 who wish to allocate their own matrix memory space. Most users should 1927 need not call this routine. 1928 1929 Level: intermediate 1930 1931 .keywords: dense, matrix, LAPACK, BLAS 1932 1933 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues(), MatSeqDenseSetLDA() 1934 1935 @*/ 1936 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 1937 { 1938 PetscErrorCode ierr,(*f)(Mat,PetscScalar[]); 1939 1940 PetscFunctionBegin; 1941 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1942 if (f) { 1943 ierr = (*f)(B,data);CHKERRQ(ierr); 1944 } 1945 PetscFunctionReturn(0); 1946 } 1947 1948 EXTERN_C_BEGIN 1949 #undef __FUNCT__ 1950 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 1951 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1952 { 1953 Mat_SeqDense *b; 1954 PetscErrorCode ierr; 1955 1956 PetscFunctionBegin; 1957 B->preallocated = PETSC_TRUE; 1958 1959 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 1960 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 1961 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1962 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1963 1964 b = (Mat_SeqDense*)B->data; 1965 b->Mmax = B->rmap->n; 1966 b->Nmax = B->cmap->n; 1967 if(b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 1968 1969 if (!data) { /* petsc-allocated storage */ 1970 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 1971 ierr = PetscMalloc(b->lda*b->Nmax*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1972 ierr = PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1973 ierr = PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1974 b->user_alloc = PETSC_FALSE; 1975 } else { /* user-allocated storage */ 1976 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 1977 b->v = data; 1978 b->user_alloc = PETSC_TRUE; 1979 } 1980 B->assembled = PETSC_TRUE; 1981 PetscFunctionReturn(0); 1982 } 1983 EXTERN_C_END 1984 1985 #undef __FUNCT__ 1986 #define __FUNCT__ "MatSeqDenseSetLDA" 1987 /*@C 1988 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 1989 1990 Input parameter: 1991 + A - the matrix 1992 - lda - the leading dimension 1993 1994 Notes: 1995 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 1996 it asserts that the preallocation has a leading dimension (the LDA parameter 1997 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 1998 1999 Level: intermediate 2000 2001 .keywords: dense, matrix, LAPACK, BLAS 2002 2003 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2004 2005 @*/ 2006 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetLDA(Mat B,PetscInt lda) 2007 { 2008 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2009 2010 PetscFunctionBegin; 2011 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); 2012 b->lda = lda; 2013 b->changelda = PETSC_FALSE; 2014 b->Mmax = PetscMax(b->Mmax,lda); 2015 PetscFunctionReturn(0); 2016 } 2017 2018 /*MC 2019 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2020 2021 Options Database Keys: 2022 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2023 2024 Level: beginner 2025 2026 .seealso: MatCreateSeqDense() 2027 2028 M*/ 2029 2030 EXTERN_C_BEGIN 2031 #undef __FUNCT__ 2032 #define __FUNCT__ "MatCreate_SeqDense" 2033 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqDense(Mat B) 2034 { 2035 Mat_SeqDense *b; 2036 PetscErrorCode ierr; 2037 PetscMPIInt size; 2038 2039 PetscFunctionBegin; 2040 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 2041 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2042 2043 ierr = PetscNewLog(B,Mat_SeqDense,&b);CHKERRQ(ierr); 2044 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2045 B->mapping = 0; 2046 B->data = (void*)b; 2047 2048 b->pivots = 0; 2049 b->roworiented = PETSC_TRUE; 2050 b->v = 0; 2051 b->changelda = PETSC_FALSE; 2052 2053 2054 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 2055 "MatGetFactor_seqdense_petsc", 2056 MatGetFactor_seqdense_petsc);CHKERRQ(ierr); 2057 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 2058 "MatSeqDenseSetPreallocation_SeqDense", 2059 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2060 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 2061 "MatMatMult_SeqAIJ_SeqDense", 2062 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2063 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C", 2064 "MatMatMultSymbolic_SeqAIJ_SeqDense", 2065 MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2066 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C", 2067 "MatMatMultNumeric_SeqAIJ_SeqDense", 2068 MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2069 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2070 PetscFunctionReturn(0); 2071 } 2072 EXTERN_C_END 2073