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