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