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