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