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_KEEP_NONZERO_PATTERN: 1358 case MAT_IGNORE_OFF_PROC_ENTRIES: 1359 case MAT_USE_HASH_TABLE: 1360 case MAT_SYMMETRIC: 1361 case MAT_STRUCTURALLY_SYMMETRIC: 1362 case MAT_HERMITIAN: 1363 case MAT_SYMMETRY_ETERNAL: 1364 case MAT_IGNORE_LOWER_TRIANGULAR: 1365 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1366 break; 1367 default: 1368 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1369 } 1370 PetscFunctionReturn(0); 1371 } 1372 1373 #undef __FUNCT__ 1374 #define __FUNCT__ "MatZeroEntries_SeqDense" 1375 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1376 { 1377 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1378 PetscErrorCode ierr; 1379 PetscInt lda=l->lda,m=A->rmap->n,j; 1380 1381 PetscFunctionBegin; 1382 if (lda>m) { 1383 for (j=0; j<A->cmap->n; j++) { 1384 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1385 } 1386 } else { 1387 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1388 } 1389 PetscFunctionReturn(0); 1390 } 1391 1392 #undef __FUNCT__ 1393 #define __FUNCT__ "MatZeroRows_SeqDense" 1394 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1395 { 1396 PetscErrorCode ierr; 1397 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1398 PetscInt m = l->lda, n = A->cmap->n, i,j; 1399 PetscScalar *slot,*bb; 1400 const PetscScalar *xx; 1401 1402 PetscFunctionBegin; 1403 #if defined(PETSC_USE_DEBUG) 1404 for (i=0; i<N; i++) { 1405 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1406 if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n); 1407 } 1408 #endif 1409 1410 /* fix right hand side if needed */ 1411 if (x && b) { 1412 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1413 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1414 for (i=0; i<N; i++) { 1415 bb[rows[i]] = diag*xx[rows[i]]; 1416 } 1417 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1418 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1419 } 1420 1421 for (i=0; i<N; i++) { 1422 slot = l->v + rows[i]; 1423 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1424 } 1425 if (diag != 0.0) { 1426 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1427 for (i=0; i<N; i++) { 1428 slot = l->v + (m+1)*rows[i]; 1429 *slot = diag; 1430 } 1431 } 1432 PetscFunctionReturn(0); 1433 } 1434 1435 #undef __FUNCT__ 1436 #define __FUNCT__ "MatGetArray_SeqDense" 1437 PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[]) 1438 { 1439 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1440 1441 PetscFunctionBegin; 1442 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"); 1443 *array = mat->v; 1444 PetscFunctionReturn(0); 1445 } 1446 1447 #undef __FUNCT__ 1448 #define __FUNCT__ "MatRestoreArray_SeqDense" 1449 PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1450 { 1451 PetscFunctionBegin; 1452 *array = 0; /* user cannot accidently use the array later */ 1453 PetscFunctionReturn(0); 1454 } 1455 1456 #undef __FUNCT__ 1457 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1458 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1459 { 1460 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1461 PetscErrorCode ierr; 1462 PetscInt i,j,nrows,ncols; 1463 const PetscInt *irow,*icol; 1464 PetscScalar *av,*bv,*v = mat->v; 1465 Mat newmat; 1466 1467 PetscFunctionBegin; 1468 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1469 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1470 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1471 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1472 1473 /* Check submatrixcall */ 1474 if (scall == MAT_REUSE_MATRIX) { 1475 PetscInt n_cols,n_rows; 1476 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1477 if (n_rows != nrows || n_cols != ncols) { 1478 /* resize the result matrix to match number of requested rows/columns */ 1479 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1480 } 1481 newmat = *B; 1482 } else { 1483 /* Create and fill new matrix */ 1484 ierr = MatCreate(((PetscObject)A)->comm,&newmat);CHKERRQ(ierr); 1485 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1486 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1487 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1488 } 1489 1490 /* Now extract the data pointers and do the copy,column at a time */ 1491 bv = ((Mat_SeqDense*)newmat->data)->v; 1492 1493 for (i=0; i<ncols; i++) { 1494 av = v + mat->lda*icol[i]; 1495 for (j=0; j<nrows; j++) { 1496 *bv++ = av[irow[j]]; 1497 } 1498 } 1499 1500 /* Assemble the matrices so that the correct flags are set */ 1501 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1502 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1503 1504 /* Free work space */ 1505 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1506 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1507 *B = newmat; 1508 PetscFunctionReturn(0); 1509 } 1510 1511 #undef __FUNCT__ 1512 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1513 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1514 { 1515 PetscErrorCode ierr; 1516 PetscInt i; 1517 1518 PetscFunctionBegin; 1519 if (scall == MAT_INITIAL_MATRIX) { 1520 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1521 } 1522 1523 for (i=0; i<n; i++) { 1524 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1525 } 1526 PetscFunctionReturn(0); 1527 } 1528 1529 #undef __FUNCT__ 1530 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1531 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1532 { 1533 PetscFunctionBegin; 1534 PetscFunctionReturn(0); 1535 } 1536 1537 #undef __FUNCT__ 1538 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1539 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1540 { 1541 PetscFunctionBegin; 1542 PetscFunctionReturn(0); 1543 } 1544 1545 #undef __FUNCT__ 1546 #define __FUNCT__ "MatCopy_SeqDense" 1547 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1548 { 1549 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1550 PetscErrorCode ierr; 1551 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1552 1553 PetscFunctionBegin; 1554 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1555 if (A->ops->copy != B->ops->copy) { 1556 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1557 PetscFunctionReturn(0); 1558 } 1559 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1560 if (lda1>m || lda2>m) { 1561 for (j=0; j<n; j++) { 1562 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1563 } 1564 } else { 1565 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1566 } 1567 PetscFunctionReturn(0); 1568 } 1569 1570 #undef __FUNCT__ 1571 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1572 PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A) 1573 { 1574 PetscErrorCode ierr; 1575 1576 PetscFunctionBegin; 1577 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1578 PetscFunctionReturn(0); 1579 } 1580 1581 #undef __FUNCT__ 1582 #define __FUNCT__ "MatSetSizes_SeqDense" 1583 PetscErrorCode MatSetSizes_SeqDense(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) 1584 { 1585 PetscFunctionBegin; 1586 /* this will not be called before lda, Mmax, and Nmax have been set */ 1587 m = PetscMax(m,M); 1588 n = PetscMax(n,N); 1589 1590 /* 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); 1591 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); 1592 */ 1593 A->rmap->n = A->rmap->N = m; 1594 A->cmap->n = A->cmap->N = n; 1595 PetscFunctionReturn(0); 1596 } 1597 1598 #undef __FUNCT__ 1599 #define __FUNCT__ "MatConjugate_SeqDense" 1600 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1601 { 1602 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1603 PetscInt i,nz = A->rmap->n*A->cmap->n; 1604 PetscScalar *aa = a->v; 1605 1606 PetscFunctionBegin; 1607 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1608 PetscFunctionReturn(0); 1609 } 1610 1611 #undef __FUNCT__ 1612 #define __FUNCT__ "MatRealPart_SeqDense" 1613 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1614 { 1615 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1616 PetscInt i,nz = A->rmap->n*A->cmap->n; 1617 PetscScalar *aa = a->v; 1618 1619 PetscFunctionBegin; 1620 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1621 PetscFunctionReturn(0); 1622 } 1623 1624 #undef __FUNCT__ 1625 #define __FUNCT__ "MatImaginaryPart_SeqDense" 1626 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 1627 { 1628 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1629 PetscInt i,nz = A->rmap->n*A->cmap->n; 1630 PetscScalar *aa = a->v; 1631 1632 PetscFunctionBegin; 1633 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1634 PetscFunctionReturn(0); 1635 } 1636 1637 /* ----------------------------------------------------------------*/ 1638 #undef __FUNCT__ 1639 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1640 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1641 { 1642 PetscErrorCode ierr; 1643 1644 PetscFunctionBegin; 1645 if (scall == MAT_INITIAL_MATRIX){ 1646 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1647 } 1648 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1649 PetscFunctionReturn(0); 1650 } 1651 1652 #undef __FUNCT__ 1653 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1654 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1655 { 1656 PetscErrorCode ierr; 1657 PetscInt m=A->rmap->n,n=B->cmap->n; 1658 Mat Cmat; 1659 1660 PetscFunctionBegin; 1661 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); 1662 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1663 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1664 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1665 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1666 Cmat->assembled = PETSC_TRUE; 1667 Cmat->ops->matmult = MatMatMult_SeqDense_SeqDense; 1668 *C = Cmat; 1669 PetscFunctionReturn(0); 1670 } 1671 1672 #undef __FUNCT__ 1673 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1674 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1675 { 1676 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1677 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1678 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1679 PetscBLASInt m,n,k; 1680 PetscScalar _DOne=1.0,_DZero=0.0; 1681 1682 PetscFunctionBegin; 1683 m = PetscBLASIntCast(A->rmap->n); 1684 n = PetscBLASIntCast(B->cmap->n); 1685 k = PetscBLASIntCast(A->cmap->n); 1686 BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1687 PetscFunctionReturn(0); 1688 } 1689 1690 #undef __FUNCT__ 1691 #define __FUNCT__ "MatTransposeMatMult_SeqDense_SeqDense" 1692 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1693 { 1694 PetscErrorCode ierr; 1695 1696 PetscFunctionBegin; 1697 if (scall == MAT_INITIAL_MATRIX){ 1698 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1699 } 1700 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1701 PetscFunctionReturn(0); 1702 } 1703 1704 #undef __FUNCT__ 1705 #define __FUNCT__ "MatTransposeMatMultSymbolic_SeqDense_SeqDense" 1706 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1707 { 1708 PetscErrorCode ierr; 1709 PetscInt m=A->cmap->n,n=B->cmap->n; 1710 Mat Cmat; 1711 1712 PetscFunctionBegin; 1713 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); 1714 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1715 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1716 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1717 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1718 Cmat->assembled = PETSC_TRUE; 1719 *C = Cmat; 1720 PetscFunctionReturn(0); 1721 } 1722 1723 #undef __FUNCT__ 1724 #define __FUNCT__ "MatTransposeMatMultNumeric_SeqDense_SeqDense" 1725 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1726 { 1727 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1728 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1729 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1730 PetscBLASInt m,n,k; 1731 PetscScalar _DOne=1.0,_DZero=0.0; 1732 1733 PetscFunctionBegin; 1734 m = PetscBLASIntCast(A->cmap->n); 1735 n = PetscBLASIntCast(B->cmap->n); 1736 k = PetscBLASIntCast(A->rmap->n); 1737 /* 1738 Note the m and n arguments below are the number rows and columns of A', not A! 1739 */ 1740 BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1741 PetscFunctionReturn(0); 1742 } 1743 1744 #undef __FUNCT__ 1745 #define __FUNCT__ "MatGetRowMax_SeqDense" 1746 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1747 { 1748 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1749 PetscErrorCode ierr; 1750 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1751 PetscScalar *x; 1752 MatScalar *aa = a->v; 1753 1754 PetscFunctionBegin; 1755 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1756 1757 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1758 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1759 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1760 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1761 for (i=0; i<m; i++) { 1762 x[i] = aa[i]; if (idx) idx[i] = 0; 1763 for (j=1; j<n; j++){ 1764 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1765 } 1766 } 1767 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1768 PetscFunctionReturn(0); 1769 } 1770 1771 #undef __FUNCT__ 1772 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1773 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1774 { 1775 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1776 PetscErrorCode ierr; 1777 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1778 PetscScalar *x; 1779 PetscReal atmp; 1780 MatScalar *aa = a->v; 1781 1782 PetscFunctionBegin; 1783 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1784 1785 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1786 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1787 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1788 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1789 for (i=0; i<m; i++) { 1790 x[i] = PetscAbsScalar(aa[i]); 1791 for (j=1; j<n; j++){ 1792 atmp = PetscAbsScalar(aa[i+m*j]); 1793 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1794 } 1795 } 1796 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1797 PetscFunctionReturn(0); 1798 } 1799 1800 #undef __FUNCT__ 1801 #define __FUNCT__ "MatGetRowMin_SeqDense" 1802 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1803 { 1804 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1805 PetscErrorCode ierr; 1806 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1807 PetscScalar *x; 1808 MatScalar *aa = a->v; 1809 1810 PetscFunctionBegin; 1811 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1812 1813 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1814 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1815 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1816 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1817 for (i=0; i<m; i++) { 1818 x[i] = aa[i]; if (idx) idx[i] = 0; 1819 for (j=1; j<n; j++){ 1820 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1821 } 1822 } 1823 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1824 PetscFunctionReturn(0); 1825 } 1826 1827 #undef __FUNCT__ 1828 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1829 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1830 { 1831 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1832 PetscErrorCode ierr; 1833 PetscScalar *x; 1834 1835 PetscFunctionBegin; 1836 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1837 1838 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1839 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1840 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1841 PetscFunctionReturn(0); 1842 } 1843 1844 1845 #undef __FUNCT__ 1846 #define __FUNCT__ "MatGetColumnNorms_SeqDense" 1847 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 1848 { 1849 PetscErrorCode ierr; 1850 PetscInt i,j,m,n; 1851 PetscScalar *a; 1852 1853 PetscFunctionBegin; 1854 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 1855 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 1856 ierr = MatGetArray(A,&a);CHKERRQ(ierr); 1857 if (type == NORM_2) { 1858 for (i=0; i<n; i++ ){ 1859 for (j=0; j<m; j++) { 1860 norms[i] += PetscAbsScalar(a[j]*a[j]); 1861 } 1862 a += m; 1863 } 1864 } else if (type == NORM_1) { 1865 for (i=0; i<n; i++ ){ 1866 for (j=0; j<m; j++) { 1867 norms[i] += PetscAbsScalar(a[j]); 1868 } 1869 a += m; 1870 } 1871 } else if (type == NORM_INFINITY) { 1872 for (i=0; i<n; i++ ){ 1873 for (j=0; j<m; j++) { 1874 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 1875 } 1876 a += m; 1877 } 1878 } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Unknown NormType"); 1879 if (type == NORM_2) { 1880 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 1881 } 1882 PetscFunctionReturn(0); 1883 } 1884 1885 /* -------------------------------------------------------------------*/ 1886 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1887 MatGetRow_SeqDense, 1888 MatRestoreRow_SeqDense, 1889 MatMult_SeqDense, 1890 /* 4*/ MatMultAdd_SeqDense, 1891 MatMultTranspose_SeqDense, 1892 MatMultTransposeAdd_SeqDense, 1893 0, 1894 0, 1895 0, 1896 /*10*/ 0, 1897 MatLUFactor_SeqDense, 1898 MatCholeskyFactor_SeqDense, 1899 MatSOR_SeqDense, 1900 MatTranspose_SeqDense, 1901 /*15*/ MatGetInfo_SeqDense, 1902 MatEqual_SeqDense, 1903 MatGetDiagonal_SeqDense, 1904 MatDiagonalScale_SeqDense, 1905 MatNorm_SeqDense, 1906 /*20*/ MatAssemblyBegin_SeqDense, 1907 MatAssemblyEnd_SeqDense, 1908 MatSetOption_SeqDense, 1909 MatZeroEntries_SeqDense, 1910 /*24*/ MatZeroRows_SeqDense, 1911 0, 1912 0, 1913 0, 1914 0, 1915 /*29*/ MatSetUpPreallocation_SeqDense, 1916 0, 1917 0, 1918 MatGetArray_SeqDense, 1919 MatRestoreArray_SeqDense, 1920 /*34*/ MatDuplicate_SeqDense, 1921 0, 1922 0, 1923 0, 1924 0, 1925 /*39*/ MatAXPY_SeqDense, 1926 MatGetSubMatrices_SeqDense, 1927 0, 1928 MatGetValues_SeqDense, 1929 MatCopy_SeqDense, 1930 /*44*/ MatGetRowMax_SeqDense, 1931 MatScale_SeqDense, 1932 0, 1933 0, 1934 0, 1935 /*49*/ 0, 1936 0, 1937 0, 1938 0, 1939 0, 1940 /*54*/ 0, 1941 0, 1942 0, 1943 0, 1944 0, 1945 /*59*/ 0, 1946 MatDestroy_SeqDense, 1947 MatView_SeqDense, 1948 0, 1949 0, 1950 /*64*/ 0, 1951 0, 1952 0, 1953 0, 1954 0, 1955 /*69*/ MatGetRowMaxAbs_SeqDense, 1956 0, 1957 0, 1958 0, 1959 0, 1960 /*74*/ 0, 1961 0, 1962 0, 1963 0, 1964 0, 1965 /*79*/ 0, 1966 0, 1967 0, 1968 0, 1969 /*83*/ MatLoad_SeqDense, 1970 0, 1971 MatIsHermitian_SeqDense, 1972 0, 1973 0, 1974 0, 1975 /*89*/ MatMatMult_SeqDense_SeqDense, 1976 MatMatMultSymbolic_SeqDense_SeqDense, 1977 MatMatMultNumeric_SeqDense_SeqDense, 1978 0, 1979 0, 1980 /*94*/ 0, 1981 0, 1982 0, 1983 0, 1984 0, 1985 /*99*/ 0, 1986 0, 1987 0, 1988 MatConjugate_SeqDense, 1989 MatSetSizes_SeqDense, 1990 /*104*/0, 1991 MatRealPart_SeqDense, 1992 MatImaginaryPart_SeqDense, 1993 0, 1994 0, 1995 /*109*/MatMatSolve_SeqDense, 1996 0, 1997 MatGetRowMin_SeqDense, 1998 MatGetColumnVector_SeqDense, 1999 0, 2000 /*114*/0, 2001 0, 2002 0, 2003 0, 2004 0, 2005 /*119*/0, 2006 0, 2007 0, 2008 0, 2009 0, 2010 /*124*/0, 2011 MatGetColumnNorms_SeqDense, 2012 0, 2013 0, 2014 0, 2015 /*129*/0, 2016 MatTransposeMatMult_SeqDense_SeqDense, 2017 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2018 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2019 }; 2020 2021 #undef __FUNCT__ 2022 #define __FUNCT__ "MatCreateSeqDense" 2023 /*@C 2024 MatCreateSeqDense - Creates a sequential dense matrix that 2025 is stored in column major order (the usual Fortran 77 manner). Many 2026 of the matrix operations use the BLAS and LAPACK routines. 2027 2028 Collective on MPI_Comm 2029 2030 Input Parameters: 2031 + comm - MPI communicator, set to PETSC_COMM_SELF 2032 . m - number of rows 2033 . n - number of columns 2034 - data - optional location of matrix data in column major order. Set data=PETSC_NULL for PETSc 2035 to control all matrix memory allocation. 2036 2037 Output Parameter: 2038 . A - the matrix 2039 2040 Notes: 2041 The data input variable is intended primarily for Fortran programmers 2042 who wish to allocate their own matrix memory space. Most users should 2043 set data=PETSC_NULL. 2044 2045 Level: intermediate 2046 2047 .keywords: dense, matrix, LAPACK, BLAS 2048 2049 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 2050 @*/ 2051 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2052 { 2053 PetscErrorCode ierr; 2054 2055 PetscFunctionBegin; 2056 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2057 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2058 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2059 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2060 PetscFunctionReturn(0); 2061 } 2062 2063 #undef __FUNCT__ 2064 #define __FUNCT__ "MatSeqDenseSetPreallocation" 2065 /*@C 2066 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2067 2068 Collective on MPI_Comm 2069 2070 Input Parameters: 2071 + A - the matrix 2072 - data - the array (or PETSC_NULL) 2073 2074 Notes: 2075 The data input variable is intended primarily for Fortran programmers 2076 who wish to allocate their own matrix memory space. Most users should 2077 need not call this routine. 2078 2079 Level: intermediate 2080 2081 .keywords: dense, matrix, LAPACK, BLAS 2082 2083 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues(), MatSeqDenseSetLDA() 2084 2085 @*/ 2086 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2087 { 2088 PetscErrorCode ierr; 2089 2090 PetscFunctionBegin; 2091 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2092 PetscFunctionReturn(0); 2093 } 2094 2095 EXTERN_C_BEGIN 2096 #undef __FUNCT__ 2097 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 2098 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2099 { 2100 Mat_SeqDense *b; 2101 PetscErrorCode ierr; 2102 2103 PetscFunctionBegin; 2104 B->preallocated = PETSC_TRUE; 2105 2106 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 2107 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 2108 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2109 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2110 2111 b = (Mat_SeqDense*)B->data; 2112 b->Mmax = B->rmap->n; 2113 b->Nmax = B->cmap->n; 2114 if(b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2115 2116 if (!data) { /* petsc-allocated storage */ 2117 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2118 ierr = PetscMalloc(b->lda*b->Nmax*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 2119 ierr = PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2120 ierr = PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2121 b->user_alloc = PETSC_FALSE; 2122 } else { /* user-allocated storage */ 2123 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2124 b->v = data; 2125 b->user_alloc = PETSC_TRUE; 2126 } 2127 B->assembled = PETSC_TRUE; 2128 PetscFunctionReturn(0); 2129 } 2130 EXTERN_C_END 2131 2132 #undef __FUNCT__ 2133 #define __FUNCT__ "MatSeqDenseSetLDA" 2134 /*@C 2135 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2136 2137 Input parameter: 2138 + A - the matrix 2139 - lda - the leading dimension 2140 2141 Notes: 2142 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2143 it asserts that the preallocation has a leading dimension (the LDA parameter 2144 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2145 2146 Level: intermediate 2147 2148 .keywords: dense, matrix, LAPACK, BLAS 2149 2150 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2151 2152 @*/ 2153 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2154 { 2155 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2156 2157 PetscFunctionBegin; 2158 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); 2159 b->lda = lda; 2160 b->changelda = PETSC_FALSE; 2161 b->Mmax = PetscMax(b->Mmax,lda); 2162 PetscFunctionReturn(0); 2163 } 2164 2165 /*MC 2166 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2167 2168 Options Database Keys: 2169 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2170 2171 Level: beginner 2172 2173 .seealso: MatCreateSeqDense() 2174 2175 M*/ 2176 2177 EXTERN_C_BEGIN 2178 #undef __FUNCT__ 2179 #define __FUNCT__ "MatCreate_SeqDense" 2180 PetscErrorCode MatCreate_SeqDense(Mat B) 2181 { 2182 Mat_SeqDense *b; 2183 PetscErrorCode ierr; 2184 PetscMPIInt size; 2185 2186 PetscFunctionBegin; 2187 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 2188 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2189 2190 ierr = PetscNewLog(B,Mat_SeqDense,&b);CHKERRQ(ierr); 2191 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2192 B->data = (void*)b; 2193 2194 b->pivots = 0; 2195 b->roworiented = PETSC_TRUE; 2196 b->v = 0; 2197 b->changelda = PETSC_FALSE; 2198 2199 2200 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 2201 "MatGetFactor_seqdense_petsc", 2202 MatGetFactor_seqdense_petsc);CHKERRQ(ierr); 2203 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 2204 "MatSeqDenseSetPreallocation_SeqDense", 2205 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2206 2207 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 2208 "MatMatMult_SeqAIJ_SeqDense", 2209 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2210 2211 2212 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C", 2213 "MatMatMultSymbolic_SeqAIJ_SeqDense", 2214 MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2215 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C", 2216 "MatMatMultNumeric_SeqAIJ_SeqDense", 2217 MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2218 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2219 PetscFunctionReturn(0); 2220 } 2221 EXTERN_C_END 2222