1 /* 2 Defines the basic matrix operations for sequential dense. 3 */ 4 5 #include "src/mat/impls/dense/seq/dense.h" 6 #include "petscblaslapack.h" 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "MatAXPY_SeqDense" 10 PetscErrorCode MatAXPY_SeqDense(const PetscScalar *alpha,Mat X,Mat Y,MatStructure str) 11 { 12 Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data; 13 int j; 14 PetscBLASInt N = (PetscBLASInt)X->m*X->n,m=(PetscBLASInt)X->m,ldax = x->lda,lday=y->lda,one = 1; 15 16 PetscFunctionBegin; 17 if (X->m != Y->m || X->n != Y->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 18 if (ldax>m || lday>m) { 19 for (j=0; j<X->n; j++) { 20 BLaxpy_(&m,(PetscScalar*)alpha,x->v+j*ldax,&one,y->v+j*lday,&one); 21 } 22 } else { 23 BLaxpy_(&N,(PetscScalar*)alpha,x->v,&one,y->v,&one); 24 } 25 PetscLogFlops(2*N-1); 26 PetscFunctionReturn(0); 27 } 28 29 #undef __FUNCT__ 30 #define __FUNCT__ "MatGetInfo_SeqDense" 31 PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 32 { 33 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 34 int i,N = A->m*A->n,count = 0; 35 PetscScalar *v = a->v; 36 37 PetscFunctionBegin; 38 for (i=0; i<N; i++) {if (*v != 0.0) count++; v++;} 39 40 info->rows_global = (double)A->m; 41 info->columns_global = (double)A->n; 42 info->rows_local = (double)A->m; 43 info->columns_local = (double)A->n; 44 info->block_size = 1.0; 45 info->nz_allocated = (double)N; 46 info->nz_used = (double)count; 47 info->nz_unneeded = (double)(N-count); 48 info->assemblies = (double)A->num_ass; 49 info->mallocs = 0; 50 info->memory = A->mem; 51 info->fill_ratio_given = 0; 52 info->fill_ratio_needed = 0; 53 info->factor_mallocs = 0; 54 55 PetscFunctionReturn(0); 56 } 57 58 #undef __FUNCT__ 59 #define __FUNCT__ "MatScale_SeqDense" 60 PetscErrorCode MatScale_SeqDense(const PetscScalar *alpha,Mat A) 61 { 62 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 63 PetscBLASInt one = 1,lda = a->lda,j,nz; 64 65 PetscFunctionBegin; 66 if (lda>A->m) { 67 nz = (PetscBLASInt)A->m; 68 for (j=0; j<A->n; j++) { 69 BLscal_(&nz,(PetscScalar*)alpha,a->v+j*lda,&one); 70 } 71 } else { 72 nz = (PetscBLASInt)A->m*A->n; 73 BLscal_(&nz,(PetscScalar*)alpha,a->v,&one); 74 } 75 PetscLogFlops(nz); 76 PetscFunctionReturn(0); 77 } 78 79 /* ---------------------------------------------------------------*/ 80 /* COMMENT: I have chosen to hide row permutation in the pivots, 81 rather than put it in the Mat->row slot.*/ 82 #undef __FUNCT__ 83 #define __FUNCT__ "MatLUFactor_SeqDense" 84 PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,MatFactorInfo *minfo) 85 { 86 #if defined(PETSC_MISSING_LAPACK_GETRF) 87 PetscFunctionBegin; 88 SETERRQ(PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable."); 89 #else 90 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 91 PetscErrorCode ierr; 92 PetscBLASInt n = (PetscBLASInt)A->n,m = (PetscBLASInt)A->m,info; 93 94 PetscFunctionBegin; 95 if (!mat->pivots) { 96 ierr = PetscMalloc((A->m+1)*sizeof(PetscBLASInt),&mat->pivots);CHKERRQ(ierr); 97 PetscLogObjectMemory(A,A->m*sizeof(PetscBLASInt)); 98 } 99 A->factor = FACTOR_LU; 100 if (!A->m || !A->n) PetscFunctionReturn(0); 101 LAgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info); 102 if (info<0) SETERRQ(PETSC_ERR_LIB,"Bad argument to LU factorization"); 103 if (info>0) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); 104 PetscLogFlops((2*A->n*A->n*A->n)/3); 105 #endif 106 PetscFunctionReturn(0); 107 } 108 109 #undef __FUNCT__ 110 #define __FUNCT__ "MatDuplicate_SeqDense" 111 PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 112 { 113 Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l; 114 PetscErrorCode ierr; 115 int lda = (int)mat->lda,j,m; 116 Mat newi; 117 118 PetscFunctionBegin; 119 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&newi);CHKERRQ(ierr); 120 ierr = MatSetType(newi,A->type_name);CHKERRQ(ierr); 121 ierr = MatSeqDenseSetPreallocation(newi,PETSC_NULL);CHKERRQ(ierr); 122 if (cpvalues == MAT_COPY_VALUES) { 123 l = (Mat_SeqDense*)newi->data; 124 if (lda>A->m) { 125 m = A->m; 126 for (j=0; j<A->n; j++) { 127 ierr = PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 128 } 129 } else { 130 ierr = PetscMemcpy(l->v,mat->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 131 } 132 } 133 newi->assembled = PETSC_TRUE; 134 *newmat = newi; 135 PetscFunctionReturn(0); 136 } 137 138 #undef __FUNCT__ 139 #define __FUNCT__ "MatLUFactorSymbolic_SeqDense" 140 PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat A,IS row,IS col,MatFactorInfo *info,Mat *fact) 141 { 142 PetscErrorCode ierr; 143 144 PetscFunctionBegin; 145 ierr = MatDuplicate_SeqDense(A,MAT_DO_NOT_COPY_VALUES,fact);CHKERRQ(ierr); 146 PetscFunctionReturn(0); 147 } 148 149 #undef __FUNCT__ 150 #define __FUNCT__ "MatLUFactorNumeric_SeqDense" 151 PetscErrorCode MatLUFactorNumeric_SeqDense(Mat A,Mat *fact) 152 { 153 Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l = (Mat_SeqDense*)(*fact)->data; 154 PetscErrorCode ierr; 155 int lda1=mat->lda,lda2=l->lda, m=A->m,n=A->n, j; 156 MatFactorInfo info; 157 158 PetscFunctionBegin; 159 /* copy the numerical values */ 160 if (lda1>m || lda2>m ) { 161 for (j=0; j<n; j++) { 162 ierr = PetscMemcpy(l->v+j*lda2,mat->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 163 } 164 } else { 165 ierr = PetscMemcpy(l->v,mat->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 166 } 167 (*fact)->factor = 0; 168 ierr = MatLUFactor(*fact,0,0,&info);CHKERRQ(ierr); 169 PetscFunctionReturn(0); 170 } 171 172 #undef __FUNCT__ 173 #define __FUNCT__ "MatCholeskyFactorSymbolic_SeqDense" 174 PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat A,IS row,MatFactorInfo *info,Mat *fact) 175 { 176 PetscErrorCode ierr; 177 178 PetscFunctionBegin; 179 ierr = MatConvert(A,MATSAME,fact);CHKERRQ(ierr); 180 PetscFunctionReturn(0); 181 } 182 183 #undef __FUNCT__ 184 #define __FUNCT__ "MatCholeskyFactor_SeqDense" 185 PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,MatFactorInfo *factinfo) 186 { 187 #if defined(PETSC_MISSING_LAPACK_POTRF) 188 PetscFunctionBegin; 189 SETERRQ(PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable."); 190 #else 191 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 192 PetscErrorCode ierr; 193 PetscBLASInt n = (PetscBLASInt)A->n,info; 194 195 PetscFunctionBegin; 196 if (mat->pivots) { 197 ierr = PetscFree(mat->pivots);CHKERRQ(ierr); 198 PetscLogObjectMemory(A,-A->m*sizeof(int)); 199 mat->pivots = 0; 200 } 201 if (!A->m || !A->n) PetscFunctionReturn(0); 202 LApotrf_("L",&n,mat->v,&mat->lda,&info); 203 if (info) SETERRQ1(PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1); 204 A->factor = FACTOR_CHOLESKY; 205 PetscLogFlops((A->n*A->n*A->n)/3); 206 #endif 207 PetscFunctionReturn(0); 208 } 209 210 #undef __FUNCT__ 211 #define __FUNCT__ "MatCholeskyFactorNumeric_SeqDense" 212 PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat A,Mat *fact) 213 { 214 PetscErrorCode ierr; 215 MatFactorInfo info; 216 217 PetscFunctionBegin; 218 info.fill = 1.0; 219 ierr = MatCholeskyFactor_SeqDense(*fact,0,&info);CHKERRQ(ierr); 220 PetscFunctionReturn(0); 221 } 222 223 #undef __FUNCT__ 224 #define __FUNCT__ "MatSolve_SeqDense" 225 PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 226 { 227 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 228 PetscErrorCode ierr; 229 PetscBLASInt m = (PetscBLASInt)A->m, one = 1,info; 230 PetscScalar *x,*y; 231 232 PetscFunctionBegin; 233 if (!A->m || !A->n) PetscFunctionReturn(0); 234 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 235 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 236 ierr = PetscMemcpy(y,x,A->m*sizeof(PetscScalar));CHKERRQ(ierr); 237 if (A->factor == FACTOR_LU) { 238 #if defined(PETSC_MISSING_LAPACK_GETRS) 239 SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 240 #else 241 LAgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 242 if (info) SETERRQ(PETSC_ERR_LIB,"GETRS - Bad solve"); 243 #endif 244 } else if (A->factor == FACTOR_CHOLESKY){ 245 #if defined(PETSC_MISSING_LAPACK_POTRS) 246 SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 247 #else 248 LApotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 249 if (info) SETERRQ(PETSC_ERR_LIB,"POTRS Bad solve"); 250 #endif 251 } 252 else SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 253 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 254 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 255 PetscLogFlops(2*A->n*A->n - A->n); 256 PetscFunctionReturn(0); 257 } 258 259 #undef __FUNCT__ 260 #define __FUNCT__ "MatSolveTranspose_SeqDense" 261 PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy) 262 { 263 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 264 PetscErrorCode ierr; 265 PetscBLASInt m = (PetscBLASInt) A->m,one = 1,info; 266 PetscScalar *x,*y; 267 268 PetscFunctionBegin; 269 if (!A->m || !A->n) PetscFunctionReturn(0); 270 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 271 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 272 ierr = PetscMemcpy(y,x,A->m*sizeof(PetscScalar));CHKERRQ(ierr); 273 /* assume if pivots exist then use LU; else Cholesky */ 274 if (mat->pivots) { 275 #if defined(PETSC_MISSING_LAPACK_GETRS) 276 SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 277 #else 278 LAgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 279 if (info) SETERRQ(PETSC_ERR_LIB,"POTRS - Bad solve"); 280 #endif 281 } else { 282 #if defined(PETSC_MISSING_LAPACK_POTRS) 283 SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 284 #else 285 LApotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 286 if (info) SETERRQ(PETSC_ERR_LIB,"POTRS - Bad solve"); 287 #endif 288 } 289 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 290 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 291 PetscLogFlops(2*A->n*A->n - A->n); 292 PetscFunctionReturn(0); 293 } 294 295 #undef __FUNCT__ 296 #define __FUNCT__ "MatSolveAdd_SeqDense" 297 PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 298 { 299 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 300 PetscErrorCode ierr; 301 PetscBLASInt m = (PetscBLASInt)A->m,one = 1,info; 302 PetscScalar *x,*y,sone = 1.0; 303 Vec tmp = 0; 304 305 PetscFunctionBegin; 306 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 307 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 308 if (!A->m || !A->n) PetscFunctionReturn(0); 309 if (yy == zz) { 310 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 311 PetscLogObjectParent(A,tmp); 312 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 313 } 314 ierr = PetscMemcpy(y,x,A->m*sizeof(PetscScalar));CHKERRQ(ierr); 315 /* assume if pivots exist then use LU; else Cholesky */ 316 if (mat->pivots) { 317 #if defined(PETSC_MISSING_LAPACK_GETRS) 318 SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 319 #else 320 LAgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 321 if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve"); 322 #endif 323 } else { 324 #if defined(PETSC_MISSING_LAPACK_POTRS) 325 SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 326 #else 327 LApotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 328 if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve"); 329 #endif 330 } 331 if (tmp) {ierr = VecAXPY(&sone,tmp,yy);CHKERRQ(ierr); ierr = VecDestroy(tmp);CHKERRQ(ierr);} 332 else {ierr = VecAXPY(&sone,zz,yy);CHKERRQ(ierr);} 333 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 334 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 335 PetscLogFlops(2*A->n*A->n); 336 PetscFunctionReturn(0); 337 } 338 339 #undef __FUNCT__ 340 #define __FUNCT__ "MatSolveTransposeAdd_SeqDense" 341 PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 342 { 343 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 344 PetscErrorCode ierr; 345 PetscBLASInt m = (PetscBLASInt)A->m,one = 1,info; 346 PetscScalar *x,*y,sone = 1.0; 347 Vec tmp; 348 349 PetscFunctionBegin; 350 if (!A->m || !A->n) PetscFunctionReturn(0); 351 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 352 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 353 if (yy == zz) { 354 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 355 PetscLogObjectParent(A,tmp); 356 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 357 } 358 ierr = PetscMemcpy(y,x,A->m*sizeof(PetscScalar));CHKERRQ(ierr); 359 /* assume if pivots exist then use LU; else Cholesky */ 360 if (mat->pivots) { 361 #if defined(PETSC_MISSING_LAPACK_GETRS) 362 SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 363 #else 364 LAgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 365 if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve"); 366 #endif 367 } else { 368 #if defined(PETSC_MISSING_LAPACK_POTRS) 369 SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 370 #else 371 LApotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 372 if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve"); 373 #endif 374 } 375 if (tmp) { 376 ierr = VecAXPY(&sone,tmp,yy);CHKERRQ(ierr); 377 ierr = VecDestroy(tmp);CHKERRQ(ierr); 378 } else { 379 ierr = VecAXPY(&sone,zz,yy);CHKERRQ(ierr); 380 } 381 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 382 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 383 PetscLogFlops(2*A->n*A->n); 384 PetscFunctionReturn(0); 385 } 386 /* ------------------------------------------------------------------*/ 387 #undef __FUNCT__ 388 #define __FUNCT__ "MatRelax_SeqDense" 389 PetscErrorCode MatRelax_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,int its,int lits,Vec xx) 390 { 391 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 392 PetscScalar *x,*b,*v = mat->v,zero = 0.0,xt; 393 PetscErrorCode ierr; 394 int m = A->m,i; 395 #if !defined(PETSC_USE_COMPLEX) 396 PetscBLASInt bm = (PetscBLASInt)m, o = 1; 397 #endif 398 399 PetscFunctionBegin; 400 if (flag & SOR_ZERO_INITIAL_GUESS) { 401 /* this is a hack fix, should have another version without the second BLdot */ 402 ierr = VecSet(&zero,xx);CHKERRQ(ierr); 403 } 404 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 405 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 406 its = its*lits; 407 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 408 while (its--) { 409 if (flag & SOR_FORWARD_SWEEP){ 410 for (i=0; i<m; i++) { 411 #if defined(PETSC_USE_COMPLEX) 412 /* cannot use BLAS dot for complex because compiler/linker is 413 not happy about returning a double complex */ 414 int _i; 415 PetscScalar sum = b[i]; 416 for (_i=0; _i<m; _i++) { 417 sum -= PetscConj(v[i+_i*m])*x[_i]; 418 } 419 xt = sum; 420 #else 421 xt = b[i] - BLdot_(&bm,v+i,&bm,x,&o); 422 #endif 423 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 424 } 425 } 426 if (flag & SOR_BACKWARD_SWEEP) { 427 for (i=m-1; i>=0; i--) { 428 #if defined(PETSC_USE_COMPLEX) 429 /* cannot use BLAS dot for complex because compiler/linker is 430 not happy about returning a double complex */ 431 int _i; 432 PetscScalar sum = b[i]; 433 for (_i=0; _i<m; _i++) { 434 sum -= PetscConj(v[i+_i*m])*x[_i]; 435 } 436 xt = sum; 437 #else 438 xt = b[i] - BLdot_(&bm,v+i,&bm,x,&o); 439 #endif 440 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 441 } 442 } 443 } 444 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 445 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 446 PetscFunctionReturn(0); 447 } 448 449 /* -----------------------------------------------------------------*/ 450 #undef __FUNCT__ 451 #define __FUNCT__ "MatMultTranspose_SeqDense" 452 PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 453 { 454 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 455 PetscScalar *v = mat->v,*x,*y; 456 PetscErrorCode ierr; 457 PetscBLASInt m = (PetscBLASInt)A->m, n = (PetscBLASInt)A->n,_One=1; 458 PetscScalar _DOne=1.0,_DZero=0.0; 459 460 PetscFunctionBegin; 461 if (!A->m || !A->n) PetscFunctionReturn(0); 462 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 463 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 464 LAgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One); 465 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 466 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 467 PetscLogFlops(2*A->m*A->n - A->n); 468 PetscFunctionReturn(0); 469 } 470 471 #undef __FUNCT__ 472 #define __FUNCT__ "MatMult_SeqDense" 473 PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 474 { 475 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 476 PetscScalar *v = mat->v,*x,*y,_DOne=1.0,_DZero=0.0; 477 PetscErrorCode ierr; 478 PetscBLASInt m = (PetscBLASInt)A->m, n = (PetscBLASInt)A->n, _One=1; 479 480 PetscFunctionBegin; 481 if (!A->m || !A->n) PetscFunctionReturn(0); 482 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 483 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 484 LAgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One); 485 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 486 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 487 PetscLogFlops(2*A->m*A->n - A->m); 488 PetscFunctionReturn(0); 489 } 490 491 #undef __FUNCT__ 492 #define __FUNCT__ "MatMultAdd_SeqDense" 493 PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 494 { 495 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 496 PetscScalar *v = mat->v,*x,*y,_DOne=1.0; 497 PetscErrorCode ierr; 498 PetscBLASInt m = (PetscBLASInt)A->m, n = (PetscBLASInt)A->n, _One=1; 499 500 PetscFunctionBegin; 501 if (!A->m || !A->n) PetscFunctionReturn(0); 502 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 503 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 504 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 505 LAgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One); 506 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 507 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 508 PetscLogFlops(2*A->m*A->n); 509 PetscFunctionReturn(0); 510 } 511 512 #undef __FUNCT__ 513 #define __FUNCT__ "MatMultTransposeAdd_SeqDense" 514 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 515 { 516 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 517 PetscScalar *v = mat->v,*x,*y; 518 PetscErrorCode ierr; 519 PetscBLASInt m = (PetscBLASInt)A->m, n = (PetscBLASInt)A->n, _One=1; 520 PetscScalar _DOne=1.0; 521 522 PetscFunctionBegin; 523 if (!A->m || !A->n) PetscFunctionReturn(0); 524 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 525 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 526 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 527 LAgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One); 528 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 529 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 530 PetscLogFlops(2*A->m*A->n); 531 PetscFunctionReturn(0); 532 } 533 534 /* -----------------------------------------------------------------*/ 535 #undef __FUNCT__ 536 #define __FUNCT__ "MatGetRow_SeqDense" 537 PetscErrorCode MatGetRow_SeqDense(Mat A,int row,int *ncols,int **cols,PetscScalar **vals) 538 { 539 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 540 PetscScalar *v; 541 PetscErrorCode ierr; 542 int i; 543 544 PetscFunctionBegin; 545 *ncols = A->n; 546 if (cols) { 547 ierr = PetscMalloc((A->n+1)*sizeof(int),cols);CHKERRQ(ierr); 548 for (i=0; i<A->n; i++) (*cols)[i] = i; 549 } 550 if (vals) { 551 ierr = PetscMalloc((A->n+1)*sizeof(PetscScalar),vals);CHKERRQ(ierr); 552 v = mat->v + row; 553 for (i=0; i<A->n; i++) {(*vals)[i] = *v; v += mat->lda;} 554 } 555 PetscFunctionReturn(0); 556 } 557 558 #undef __FUNCT__ 559 #define __FUNCT__ "MatRestoreRow_SeqDense" 560 PetscErrorCode MatRestoreRow_SeqDense(Mat A,int row,int *ncols,int **cols,PetscScalar **vals) 561 { 562 PetscErrorCode ierr; 563 PetscFunctionBegin; 564 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 565 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 566 PetscFunctionReturn(0); 567 } 568 /* ----------------------------------------------------------------*/ 569 #undef __FUNCT__ 570 #define __FUNCT__ "MatSetValues_SeqDense" 571 PetscErrorCode MatSetValues_SeqDense(Mat A,int m,const int indexm[],int n,const int indexn[],const PetscScalar v[],InsertMode addv) 572 { 573 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 574 int i,j,idx=0; 575 576 PetscFunctionBegin; 577 if (!mat->roworiented) { 578 if (addv == INSERT_VALUES) { 579 for (j=0; j<n; j++) { 580 if (indexn[j] < 0) {idx += m; continue;} 581 #if defined(PETSC_USE_BOPT_g) 582 if (indexn[j] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->n-1); 583 #endif 584 for (i=0; i<m; i++) { 585 if (indexm[i] < 0) {idx++; continue;} 586 #if defined(PETSC_USE_BOPT_g) 587 if (indexm[i] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->m-1); 588 #endif 589 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 590 } 591 } 592 } else { 593 for (j=0; j<n; j++) { 594 if (indexn[j] < 0) {idx += m; continue;} 595 #if defined(PETSC_USE_BOPT_g) 596 if (indexn[j] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->n-1); 597 #endif 598 for (i=0; i<m; i++) { 599 if (indexm[i] < 0) {idx++; continue;} 600 #if defined(PETSC_USE_BOPT_g) 601 if (indexm[i] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->m-1); 602 #endif 603 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 604 } 605 } 606 } 607 } else { 608 if (addv == INSERT_VALUES) { 609 for (i=0; i<m; i++) { 610 if (indexm[i] < 0) { idx += n; continue;} 611 #if defined(PETSC_USE_BOPT_g) 612 if (indexm[i] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->m-1); 613 #endif 614 for (j=0; j<n; j++) { 615 if (indexn[j] < 0) { idx++; continue;} 616 #if defined(PETSC_USE_BOPT_g) 617 if (indexn[j] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->n-1); 618 #endif 619 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 620 } 621 } 622 } else { 623 for (i=0; i<m; i++) { 624 if (indexm[i] < 0) { idx += n; continue;} 625 #if defined(PETSC_USE_BOPT_g) 626 if (indexm[i] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->m-1); 627 #endif 628 for (j=0; j<n; j++) { 629 if (indexn[j] < 0) { idx++; continue;} 630 #if defined(PETSC_USE_BOPT_g) 631 if (indexn[j] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->n-1); 632 #endif 633 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 634 } 635 } 636 } 637 } 638 PetscFunctionReturn(0); 639 } 640 641 #undef __FUNCT__ 642 #define __FUNCT__ "MatGetValues_SeqDense" 643 PetscErrorCode MatGetValues_SeqDense(Mat A,int m,const int indexm[],int n,const int indexn[],PetscScalar v[]) 644 { 645 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 646 int i,j; 647 PetscScalar *vpt = v; 648 649 PetscFunctionBegin; 650 /* row-oriented output */ 651 for (i=0; i<m; i++) { 652 for (j=0; j<n; j++) { 653 *vpt++ = mat->v[indexn[j]*mat->lda + indexm[i]]; 654 } 655 } 656 PetscFunctionReturn(0); 657 } 658 659 /* -----------------------------------------------------------------*/ 660 661 #include "petscsys.h" 662 663 #undef __FUNCT__ 664 #define __FUNCT__ "MatLoad_SeqDense" 665 PetscErrorCode MatLoad_SeqDense(PetscViewer viewer,const MatType type,Mat *A) 666 { 667 Mat_SeqDense *a; 668 Mat B; 669 PetscErrorCode ierr; 670 int *scols,i,j,nz,fd,header[4],size; 671 int *rowlengths = 0,M,N,*cols; 672 PetscScalar *vals,*svals,*v,*w; 673 MPI_Comm comm = ((PetscObject)viewer)->comm; 674 675 PetscFunctionBegin; 676 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 677 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 678 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 679 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 680 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 681 M = header[1]; N = header[2]; nz = header[3]; 682 683 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 684 ierr = MatCreate(comm,M,N,M,N,A);CHKERRQ(ierr); 685 ierr = MatSetType(*A,type);CHKERRQ(ierr); 686 ierr = MatSeqDenseSetPreallocation(*A,PETSC_NULL);CHKERRQ(ierr); 687 B = *A; 688 a = (Mat_SeqDense*)B->data; 689 v = a->v; 690 /* Allocate some temp space to read in the values and then flip them 691 from row major to column major */ 692 ierr = PetscMalloc((M*N > 0 ? M*N : 1)*sizeof(PetscScalar),&w);CHKERRQ(ierr); 693 /* read in nonzero values */ 694 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 695 /* now flip the values and store them in the matrix*/ 696 for (j=0; j<N; j++) { 697 for (i=0; i<M; i++) { 698 *v++ =w[i*N+j]; 699 } 700 } 701 ierr = PetscFree(w);CHKERRQ(ierr); 702 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 703 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 704 } else { 705 /* read row lengths */ 706 ierr = PetscMalloc((M+1)*sizeof(int),&rowlengths);CHKERRQ(ierr); 707 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 708 709 /* create our matrix */ 710 ierr = MatCreate(comm,M,N,M,N,A);CHKERRQ(ierr); 711 ierr = MatSetType(*A,type);CHKERRQ(ierr); 712 ierr = MatSeqDenseSetPreallocation(*A,PETSC_NULL);CHKERRQ(ierr); 713 B = *A; 714 a = (Mat_SeqDense*)B->data; 715 v = a->v; 716 717 /* read column indices and nonzeros */ 718 ierr = PetscMalloc((nz+1)*sizeof(int),&scols);CHKERRQ(ierr); 719 cols = scols; 720 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 721 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&svals);CHKERRQ(ierr); 722 vals = svals; 723 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 724 725 /* insert into matrix */ 726 for (i=0; i<M; i++) { 727 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 728 svals += rowlengths[i]; scols += rowlengths[i]; 729 } 730 ierr = PetscFree(vals);CHKERRQ(ierr); 731 ierr = PetscFree(cols);CHKERRQ(ierr); 732 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 733 734 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 735 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 736 } 737 PetscFunctionReturn(0); 738 } 739 740 #include "petscsys.h" 741 742 #undef __FUNCT__ 743 #define __FUNCT__ "MatView_SeqDense_ASCII" 744 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 745 { 746 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 747 PetscErrorCode ierr; 748 int i,j; 749 char *name; 750 PetscScalar *v; 751 PetscViewerFormat format; 752 753 PetscFunctionBegin; 754 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 755 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 756 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 757 PetscFunctionReturn(0); /* do nothing for now */ 758 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 759 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 760 for (i=0; i<A->m; i++) { 761 v = a->v + i; 762 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 763 for (j=0; j<A->n; j++) { 764 #if defined(PETSC_USE_COMPLEX) 765 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 766 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 767 } else if (PetscRealPart(*v)) { 768 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,PetscRealPart(*v));CHKERRQ(ierr); 769 } 770 #else 771 if (*v) { 772 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,*v);CHKERRQ(ierr); 773 } 774 #endif 775 v += a->lda; 776 } 777 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 778 } 779 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 780 } else { 781 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 782 #if defined(PETSC_USE_COMPLEX) 783 PetscTruth allreal = PETSC_TRUE; 784 /* determine if matrix has all real values */ 785 v = a->v; 786 for (i=0; i<A->m*A->n; i++) { 787 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;} 788 } 789 #endif 790 if (format == PETSC_VIEWER_ASCII_MATLAB) { 791 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 792 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->m,A->n);CHKERRQ(ierr); 793 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->m,A->n);CHKERRQ(ierr); 794 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 795 } 796 797 for (i=0; i<A->m; i++) { 798 v = a->v + i; 799 for (j=0; j<A->n; j++) { 800 #if defined(PETSC_USE_COMPLEX) 801 if (allreal) { 802 ierr = PetscViewerASCIIPrintf(viewer,"%6.4e ",PetscRealPart(*v));CHKERRQ(ierr); 803 } else { 804 ierr = PetscViewerASCIIPrintf(viewer,"%6.4e + %6.4e i ",PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 805 } 806 #else 807 ierr = PetscViewerASCIIPrintf(viewer,"%6.4e ",*v);CHKERRQ(ierr); 808 #endif 809 v += a->lda; 810 } 811 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 812 } 813 if (format == PETSC_VIEWER_ASCII_MATLAB) { 814 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 815 } 816 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 817 } 818 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 819 PetscFunctionReturn(0); 820 } 821 822 #undef __FUNCT__ 823 #define __FUNCT__ "MatView_SeqDense_Binary" 824 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 825 { 826 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 827 PetscErrorCode ierr; 828 PetscInt ict,j,n = A->n,m = A->m,i,fd,*col_lens,nz = m*n; 829 PetscScalar *v,*anonz,*vals; 830 PetscViewerFormat format; 831 832 PetscFunctionBegin; 833 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 834 835 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 836 if (format == PETSC_VIEWER_BINARY_NATIVE) { 837 /* store the matrix as a dense matrix */ 838 ierr = PetscMalloc(4*sizeof(int),&col_lens);CHKERRQ(ierr); 839 col_lens[0] = MAT_FILE_COOKIE; 840 col_lens[1] = m; 841 col_lens[2] = n; 842 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 843 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 844 ierr = PetscFree(col_lens);CHKERRQ(ierr); 845 846 /* write out matrix, by rows */ 847 ierr = PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 848 v = a->v; 849 for (i=0; i<m; i++) { 850 for (j=0; j<n; j++) { 851 vals[i + j*m] = *v++; 852 } 853 } 854 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 855 ierr = PetscFree(vals);CHKERRQ(ierr); 856 } else { 857 ierr = PetscMalloc((4+nz)*sizeof(int),&col_lens);CHKERRQ(ierr); 858 col_lens[0] = MAT_FILE_COOKIE; 859 col_lens[1] = m; 860 col_lens[2] = n; 861 col_lens[3] = nz; 862 863 /* store lengths of each row and write (including header) to file */ 864 for (i=0; i<m; i++) col_lens[4+i] = n; 865 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 866 867 /* Possibly should write in smaller increments, not whole matrix at once? */ 868 /* store column indices (zero start index) */ 869 ict = 0; 870 for (i=0; i<m; i++) { 871 for (j=0; j<n; j++) col_lens[ict++] = j; 872 } 873 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 874 ierr = PetscFree(col_lens);CHKERRQ(ierr); 875 876 /* store nonzero values */ 877 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);CHKERRQ(ierr); 878 ict = 0; 879 for (i=0; i<m; i++) { 880 v = a->v + i; 881 for (j=0; j<n; j++) { 882 anonz[ict++] = *v; v += a->lda; 883 } 884 } 885 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 886 ierr = PetscFree(anonz);CHKERRQ(ierr); 887 } 888 PetscFunctionReturn(0); 889 } 890 891 #undef __FUNCT__ 892 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 893 PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 894 { 895 Mat A = (Mat) Aa; 896 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 897 PetscErrorCode ierr; 898 int m = A->m,n = A->n,color,i,j; 899 PetscScalar *v = a->v; 900 PetscViewer viewer; 901 PetscDraw popup; 902 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 903 PetscViewerFormat format; 904 905 PetscFunctionBegin; 906 907 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 908 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 909 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 910 911 /* Loop over matrix elements drawing boxes */ 912 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 913 /* Blue for negative and Red for positive */ 914 color = PETSC_DRAW_BLUE; 915 for(j = 0; j < n; j++) { 916 x_l = j; 917 x_r = x_l + 1.0; 918 for(i = 0; i < m; i++) { 919 y_l = m - i - 1.0; 920 y_r = y_l + 1.0; 921 #if defined(PETSC_USE_COMPLEX) 922 if (PetscRealPart(v[j*m+i]) > 0.) { 923 color = PETSC_DRAW_RED; 924 } else if (PetscRealPart(v[j*m+i]) < 0.) { 925 color = PETSC_DRAW_BLUE; 926 } else { 927 continue; 928 } 929 #else 930 if (v[j*m+i] > 0.) { 931 color = PETSC_DRAW_RED; 932 } else if (v[j*m+i] < 0.) { 933 color = PETSC_DRAW_BLUE; 934 } else { 935 continue; 936 } 937 #endif 938 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 939 } 940 } 941 } else { 942 /* use contour shading to indicate magnitude of values */ 943 /* first determine max of all nonzero values */ 944 for(i = 0; i < m*n; i++) { 945 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 946 } 947 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 948 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 949 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 950 for(j = 0; j < n; j++) { 951 x_l = j; 952 x_r = x_l + 1.0; 953 for(i = 0; i < m; i++) { 954 y_l = m - i - 1.0; 955 y_r = y_l + 1.0; 956 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 957 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 958 } 959 } 960 } 961 PetscFunctionReturn(0); 962 } 963 964 #undef __FUNCT__ 965 #define __FUNCT__ "MatView_SeqDense_Draw" 966 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 967 { 968 PetscDraw draw; 969 PetscTruth isnull; 970 PetscReal xr,yr,xl,yl,h,w; 971 PetscErrorCode ierr; 972 973 PetscFunctionBegin; 974 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 975 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 976 if (isnull == PETSC_TRUE) PetscFunctionReturn(0); 977 978 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 979 xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0; 980 xr += w; yr += h; xl = -w; yl = -h; 981 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 982 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 983 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 984 PetscFunctionReturn(0); 985 } 986 987 #undef __FUNCT__ 988 #define __FUNCT__ "MatView_SeqDense" 989 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 990 { 991 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 992 PetscErrorCode ierr; 993 PetscTruth issocket,iascii,isbinary,isdraw; 994 995 PetscFunctionBegin; 996 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 997 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 998 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 999 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 1000 1001 if (issocket) { 1002 if (a->lda>A->m) SETERRQ(PETSC_ERR_SUP,"Case can not handle LDA"); 1003 ierr = PetscViewerSocketPutScalar(viewer,A->m,A->n,a->v);CHKERRQ(ierr); 1004 } else if (iascii) { 1005 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1006 } else if (isbinary) { 1007 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1008 } else if (isdraw) { 1009 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1010 } else { 1011 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name); 1012 } 1013 PetscFunctionReturn(0); 1014 } 1015 1016 #undef __FUNCT__ 1017 #define __FUNCT__ "MatDestroy_SeqDense" 1018 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1019 { 1020 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1021 PetscErrorCode ierr; 1022 1023 PetscFunctionBegin; 1024 #if defined(PETSC_USE_LOG) 1025 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->m,mat->n); 1026 #endif 1027 if (l->pivots) {ierr = PetscFree(l->pivots);CHKERRQ(ierr);} 1028 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1029 ierr = PetscFree(l);CHKERRQ(ierr); 1030 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 1031 PetscFunctionReturn(0); 1032 } 1033 1034 #undef __FUNCT__ 1035 #define __FUNCT__ "MatTranspose_SeqDense" 1036 PetscErrorCode MatTranspose_SeqDense(Mat A,Mat *matout) 1037 { 1038 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1039 PetscErrorCode ierr; 1040 int k,j,m,n,M; 1041 PetscScalar *v,tmp; 1042 1043 PetscFunctionBegin; 1044 v = mat->v; m = A->m; M = mat->lda; n = A->n; 1045 if (!matout) { /* in place transpose */ 1046 if (m != n) { 1047 SETERRQ(PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1048 } else { 1049 for (j=0; j<m; j++) { 1050 for (k=0; k<j; k++) { 1051 tmp = v[j + k*M]; 1052 v[j + k*M] = v[k + j*M]; 1053 v[k + j*M] = tmp; 1054 } 1055 } 1056 } 1057 } else { /* out-of-place transpose */ 1058 Mat tmat; 1059 Mat_SeqDense *tmatd; 1060 PetscScalar *v2; 1061 1062 ierr = MatCreate(A->comm,A->n,A->m,A->n,A->m,&tmat);CHKERRQ(ierr); 1063 ierr = MatSetType(tmat,A->type_name);CHKERRQ(ierr); 1064 ierr = MatSeqDenseSetPreallocation(tmat,PETSC_NULL);CHKERRQ(ierr); 1065 tmatd = (Mat_SeqDense*)tmat->data; 1066 v = mat->v; v2 = tmatd->v; 1067 for (j=0; j<n; j++) { 1068 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1069 } 1070 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1071 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1072 *matout = tmat; 1073 } 1074 PetscFunctionReturn(0); 1075 } 1076 1077 #undef __FUNCT__ 1078 #define __FUNCT__ "MatEqual_SeqDense" 1079 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg) 1080 { 1081 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1082 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1083 int i,j; 1084 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1085 1086 PetscFunctionBegin; 1087 if (A1->m != A2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1088 if (A1->n != A2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1089 for (i=0; i<A1->m; i++) { 1090 v1 = mat1->v+i; v2 = mat2->v+i; 1091 for (j=0; j<A1->n; j++) { 1092 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1093 v1 += mat1->lda; v2 += mat2->lda; 1094 } 1095 } 1096 *flg = PETSC_TRUE; 1097 PetscFunctionReturn(0); 1098 } 1099 1100 #undef __FUNCT__ 1101 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1102 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1103 { 1104 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1105 PetscErrorCode ierr; 1106 int i,n,len; 1107 PetscScalar *x,zero = 0.0; 1108 1109 PetscFunctionBegin; 1110 ierr = VecSet(&zero,v);CHKERRQ(ierr); 1111 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1112 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1113 len = PetscMin(A->m,A->n); 1114 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1115 for (i=0; i<len; i++) { 1116 x[i] = mat->v[i*mat->lda + i]; 1117 } 1118 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1119 PetscFunctionReturn(0); 1120 } 1121 1122 #undef __FUNCT__ 1123 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1124 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1125 { 1126 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1127 PetscScalar *l,*r,x,*v; 1128 PetscErrorCode ierr; 1129 int i,j,m = A->m,n = A->n; 1130 1131 PetscFunctionBegin; 1132 if (ll) { 1133 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1134 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1135 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1136 for (i=0; i<m; i++) { 1137 x = l[i]; 1138 v = mat->v + i; 1139 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1140 } 1141 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1142 PetscLogFlops(n*m); 1143 } 1144 if (rr) { 1145 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1146 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1147 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1148 for (i=0; i<n; i++) { 1149 x = r[i]; 1150 v = mat->v + i*m; 1151 for (j=0; j<m; j++) { (*v++) *= x;} 1152 } 1153 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1154 PetscLogFlops(n*m); 1155 } 1156 PetscFunctionReturn(0); 1157 } 1158 1159 #undef __FUNCT__ 1160 #define __FUNCT__ "MatNorm_SeqDense" 1161 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1162 { 1163 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1164 PetscScalar *v = mat->v; 1165 PetscReal sum = 0.0; 1166 int lda=mat->lda,m=A->m,i,j; 1167 1168 PetscFunctionBegin; 1169 if (type == NORM_FROBENIUS) { 1170 if (lda>m) { 1171 for (j=0; j<A->n; j++) { 1172 v = mat->v+j*lda; 1173 for (i=0; i<m; i++) { 1174 #if defined(PETSC_USE_COMPLEX) 1175 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1176 #else 1177 sum += (*v)*(*v); v++; 1178 #endif 1179 } 1180 } 1181 } else { 1182 for (i=0; i<A->n*A->m; i++) { 1183 #if defined(PETSC_USE_COMPLEX) 1184 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1185 #else 1186 sum += (*v)*(*v); v++; 1187 #endif 1188 } 1189 } 1190 *nrm = sqrt(sum); 1191 PetscLogFlops(2*A->n*A->m); 1192 } else if (type == NORM_1) { 1193 *nrm = 0.0; 1194 for (j=0; j<A->n; j++) { 1195 v = mat->v + j*mat->lda; 1196 sum = 0.0; 1197 for (i=0; i<A->m; i++) { 1198 sum += PetscAbsScalar(*v); v++; 1199 } 1200 if (sum > *nrm) *nrm = sum; 1201 } 1202 PetscLogFlops(A->n*A->m); 1203 } else if (type == NORM_INFINITY) { 1204 *nrm = 0.0; 1205 for (j=0; j<A->m; j++) { 1206 v = mat->v + j; 1207 sum = 0.0; 1208 for (i=0; i<A->n; i++) { 1209 sum += PetscAbsScalar(*v); v += mat->lda; 1210 } 1211 if (sum > *nrm) *nrm = sum; 1212 } 1213 PetscLogFlops(A->n*A->m); 1214 } else { 1215 SETERRQ(PETSC_ERR_SUP,"No two norm"); 1216 } 1217 PetscFunctionReturn(0); 1218 } 1219 1220 #undef __FUNCT__ 1221 #define __FUNCT__ "MatSetOption_SeqDense" 1222 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op) 1223 { 1224 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1225 1226 PetscFunctionBegin; 1227 switch (op) { 1228 case MAT_ROW_ORIENTED: 1229 aij->roworiented = PETSC_TRUE; 1230 break; 1231 case MAT_COLUMN_ORIENTED: 1232 aij->roworiented = PETSC_FALSE; 1233 break; 1234 case MAT_ROWS_SORTED: 1235 case MAT_ROWS_UNSORTED: 1236 case MAT_COLUMNS_SORTED: 1237 case MAT_COLUMNS_UNSORTED: 1238 case MAT_NO_NEW_NONZERO_LOCATIONS: 1239 case MAT_YES_NEW_NONZERO_LOCATIONS: 1240 case MAT_NEW_NONZERO_LOCATION_ERR: 1241 case MAT_NO_NEW_DIAGONALS: 1242 case MAT_YES_NEW_DIAGONALS: 1243 case MAT_IGNORE_OFF_PROC_ENTRIES: 1244 case MAT_USE_HASH_TABLE: 1245 PetscLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1246 break; 1247 case MAT_SYMMETRIC: 1248 case MAT_STRUCTURALLY_SYMMETRIC: 1249 case MAT_NOT_SYMMETRIC: 1250 case MAT_NOT_STRUCTURALLY_SYMMETRIC: 1251 case MAT_HERMITIAN: 1252 case MAT_NOT_HERMITIAN: 1253 case MAT_SYMMETRY_ETERNAL: 1254 case MAT_NOT_SYMMETRY_ETERNAL: 1255 break; 1256 default: 1257 SETERRQ(PETSC_ERR_SUP,"unknown option"); 1258 } 1259 PetscFunctionReturn(0); 1260 } 1261 1262 #undef __FUNCT__ 1263 #define __FUNCT__ "MatZeroEntries_SeqDense" 1264 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1265 { 1266 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1267 PetscErrorCode ierr; 1268 int lda=l->lda,m=A->m,j; 1269 1270 PetscFunctionBegin; 1271 if (lda>m) { 1272 for (j=0; j<A->n; j++) { 1273 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1274 } 1275 } else { 1276 ierr = PetscMemzero(l->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 1277 } 1278 PetscFunctionReturn(0); 1279 } 1280 1281 #undef __FUNCT__ 1282 #define __FUNCT__ "MatZeroRows_SeqDense" 1283 PetscErrorCode MatZeroRows_SeqDense(Mat A,IS is,const PetscScalar *diag) 1284 { 1285 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1286 PetscErrorCode ierr; 1287 int n = A->n,i,j,N,*rows; 1288 PetscScalar *slot; 1289 1290 PetscFunctionBegin; 1291 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1292 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1293 for (i=0; i<N; i++) { 1294 slot = l->v + rows[i]; 1295 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1296 } 1297 if (diag) { 1298 for (i=0; i<N; i++) { 1299 slot = l->v + (n+1)*rows[i]; 1300 *slot = *diag; 1301 } 1302 } 1303 ISRestoreIndices(is,&rows); 1304 PetscFunctionReturn(0); 1305 } 1306 1307 #undef __FUNCT__ 1308 #define __FUNCT__ "MatGetArray_SeqDense" 1309 PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[]) 1310 { 1311 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1312 1313 PetscFunctionBegin; 1314 *array = mat->v; 1315 PetscFunctionReturn(0); 1316 } 1317 1318 #undef __FUNCT__ 1319 #define __FUNCT__ "MatRestoreArray_SeqDense" 1320 PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1321 { 1322 PetscFunctionBegin; 1323 *array = 0; /* user cannot accidently use the array later */ 1324 PetscFunctionReturn(0); 1325 } 1326 1327 #undef __FUNCT__ 1328 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1329 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatReuse scall,Mat *B) 1330 { 1331 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1332 PetscErrorCode ierr; 1333 int i,j,m = A->m,*irow,*icol,nrows,ncols; 1334 PetscScalar *av,*bv,*v = mat->v; 1335 Mat newmat; 1336 1337 PetscFunctionBegin; 1338 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1339 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1340 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1341 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1342 1343 /* Check submatrixcall */ 1344 if (scall == MAT_REUSE_MATRIX) { 1345 int n_cols,n_rows; 1346 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1347 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1348 newmat = *B; 1349 } else { 1350 /* Create and fill new matrix */ 1351 ierr = MatCreate(A->comm,nrows,ncols,nrows,ncols,&newmat);CHKERRQ(ierr); 1352 ierr = MatSetType(newmat,A->type_name);CHKERRQ(ierr); 1353 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1354 } 1355 1356 /* Now extract the data pointers and do the copy,column at a time */ 1357 bv = ((Mat_SeqDense*)newmat->data)->v; 1358 1359 for (i=0; i<ncols; i++) { 1360 av = v + m*icol[i]; 1361 for (j=0; j<nrows; j++) { 1362 *bv++ = av[irow[j]]; 1363 } 1364 } 1365 1366 /* Assemble the matrices so that the correct flags are set */ 1367 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1368 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1369 1370 /* Free work space */ 1371 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1372 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1373 *B = newmat; 1374 PetscFunctionReturn(0); 1375 } 1376 1377 #undef __FUNCT__ 1378 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1379 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,int n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1380 { 1381 PetscErrorCode ierr; 1382 int i; 1383 1384 PetscFunctionBegin; 1385 if (scall == MAT_INITIAL_MATRIX) { 1386 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1387 } 1388 1389 for (i=0; i<n; i++) { 1390 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1391 } 1392 PetscFunctionReturn(0); 1393 } 1394 1395 #undef __FUNCT__ 1396 #define __FUNCT__ "MatCopy_SeqDense" 1397 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1398 { 1399 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1400 PetscErrorCode ierr; 1401 int lda1=a->lda,lda2=b->lda, m=A->m,n=A->n, j; 1402 1403 PetscFunctionBegin; 1404 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1405 if (A->ops->copy != B->ops->copy) { 1406 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1407 PetscFunctionReturn(0); 1408 } 1409 if (m != B->m || n != B->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1410 if (lda1>m || lda2>m) { 1411 for (j=0; j<n; j++) { 1412 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1413 } 1414 } else { 1415 ierr = PetscMemcpy(b->v,a->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 1416 } 1417 PetscFunctionReturn(0); 1418 } 1419 1420 #undef __FUNCT__ 1421 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1422 PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A) 1423 { 1424 PetscErrorCode ierr; 1425 1426 PetscFunctionBegin; 1427 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1428 PetscFunctionReturn(0); 1429 } 1430 1431 /* -------------------------------------------------------------------*/ 1432 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1433 MatGetRow_SeqDense, 1434 MatRestoreRow_SeqDense, 1435 MatMult_SeqDense, 1436 /* 4*/ MatMultAdd_SeqDense, 1437 MatMultTranspose_SeqDense, 1438 MatMultTransposeAdd_SeqDense, 1439 MatSolve_SeqDense, 1440 MatSolveAdd_SeqDense, 1441 MatSolveTranspose_SeqDense, 1442 /*10*/ MatSolveTransposeAdd_SeqDense, 1443 MatLUFactor_SeqDense, 1444 MatCholeskyFactor_SeqDense, 1445 MatRelax_SeqDense, 1446 MatTranspose_SeqDense, 1447 /*15*/ MatGetInfo_SeqDense, 1448 MatEqual_SeqDense, 1449 MatGetDiagonal_SeqDense, 1450 MatDiagonalScale_SeqDense, 1451 MatNorm_SeqDense, 1452 /*20*/ 0, 1453 0, 1454 0, 1455 MatSetOption_SeqDense, 1456 MatZeroEntries_SeqDense, 1457 /*25*/ MatZeroRows_SeqDense, 1458 MatLUFactorSymbolic_SeqDense, 1459 MatLUFactorNumeric_SeqDense, 1460 MatCholeskyFactorSymbolic_SeqDense, 1461 MatCholeskyFactorNumeric_SeqDense, 1462 /*30*/ MatSetUpPreallocation_SeqDense, 1463 0, 1464 0, 1465 MatGetArray_SeqDense, 1466 MatRestoreArray_SeqDense, 1467 /*35*/ MatDuplicate_SeqDense, 1468 0, 1469 0, 1470 0, 1471 0, 1472 /*40*/ MatAXPY_SeqDense, 1473 MatGetSubMatrices_SeqDense, 1474 0, 1475 MatGetValues_SeqDense, 1476 MatCopy_SeqDense, 1477 /*45*/ 0, 1478 MatScale_SeqDense, 1479 0, 1480 0, 1481 0, 1482 /*50*/ 0, 1483 0, 1484 0, 1485 0, 1486 0, 1487 /*55*/ 0, 1488 0, 1489 0, 1490 0, 1491 0, 1492 /*60*/ 0, 1493 MatDestroy_SeqDense, 1494 MatView_SeqDense, 1495 MatGetPetscMaps_Petsc, 1496 0, 1497 /*65*/ 0, 1498 0, 1499 0, 1500 0, 1501 0, 1502 /*70*/ 0, 1503 0, 1504 0, 1505 0, 1506 0, 1507 /*75*/ 0, 1508 0, 1509 0, 1510 0, 1511 0, 1512 /*80*/ 0, 1513 0, 1514 0, 1515 0, 1516 /*84*/ MatLoad_SeqDense, 1517 0, 1518 0, 1519 0, 1520 0, 1521 0, 1522 /*90*/ 0, 1523 0, 1524 0, 1525 0, 1526 0, 1527 /*95*/ 0, 1528 0, 1529 0, 1530 0}; 1531 1532 #undef __FUNCT__ 1533 #define __FUNCT__ "MatCreateSeqDense" 1534 /*@C 1535 MatCreateSeqDense - Creates a sequential dense matrix that 1536 is stored in column major order (the usual Fortran 77 manner). Many 1537 of the matrix operations use the BLAS and LAPACK routines. 1538 1539 Collective on MPI_Comm 1540 1541 Input Parameters: 1542 + comm - MPI communicator, set to PETSC_COMM_SELF 1543 . m - number of rows 1544 . n - number of columns 1545 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1546 to control all matrix memory allocation. 1547 1548 Output Parameter: 1549 . A - the matrix 1550 1551 Notes: 1552 The data input variable is intended primarily for Fortran programmers 1553 who wish to allocate their own matrix memory space. Most users should 1554 set data=PETSC_NULL. 1555 1556 Level: intermediate 1557 1558 .keywords: dense, matrix, LAPACK, BLAS 1559 1560 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1561 @*/ 1562 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,int m,int n,PetscScalar *data,Mat *A) 1563 { 1564 PetscErrorCode ierr; 1565 1566 PetscFunctionBegin; 1567 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1568 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1569 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1570 PetscFunctionReturn(0); 1571 } 1572 1573 #undef __FUNCT__ 1574 #define __FUNCT__ "MatSeqDensePreallocation" 1575 /*@C 1576 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1577 1578 Collective on MPI_Comm 1579 1580 Input Parameters: 1581 + A - the matrix 1582 - data - the array (or PETSC_NULL) 1583 1584 Notes: 1585 The data input variable is intended primarily for Fortran programmers 1586 who wish to allocate their own matrix memory space. Most users should 1587 set data=PETSC_NULL. 1588 1589 Level: intermediate 1590 1591 .keywords: dense, matrix, LAPACK, BLAS 1592 1593 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1594 @*/ 1595 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 1596 { 1597 PetscErrorCode ierr,(*f)(Mat,PetscScalar[]); 1598 1599 PetscFunctionBegin; 1600 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1601 if (f) { 1602 ierr = (*f)(B,data);CHKERRQ(ierr); 1603 } 1604 PetscFunctionReturn(0); 1605 } 1606 1607 EXTERN_C_BEGIN 1608 #undef __FUNCT__ 1609 #define __FUNCT__ "MatSeqDensePreallocation_SeqDense" 1610 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1611 { 1612 Mat_SeqDense *b; 1613 PetscErrorCode ierr; 1614 1615 PetscFunctionBegin; 1616 B->preallocated = PETSC_TRUE; 1617 b = (Mat_SeqDense*)B->data; 1618 if (!data) { 1619 ierr = PetscMalloc((B->m*B->n+1)*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1620 ierr = PetscMemzero(b->v,B->m*B->n*sizeof(PetscScalar));CHKERRQ(ierr); 1621 b->user_alloc = PETSC_FALSE; 1622 PetscLogObjectMemory(B,B->n*B->m*sizeof(PetscScalar)); 1623 } else { /* user-allocated storage */ 1624 b->v = data; 1625 b->user_alloc = PETSC_TRUE; 1626 } 1627 PetscFunctionReturn(0); 1628 } 1629 EXTERN_C_END 1630 1631 #undef __FUNCT__ 1632 #define __FUNCT__ "MatSeqDenseSetLDA" 1633 /*@C 1634 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 1635 1636 Input parameter: 1637 + A - the matrix 1638 - lda - the leading dimension 1639 1640 Notes: 1641 This routine is to be used in conjunction with MatSeqDenseSetPreallocation; 1642 it asserts that the preallocation has a leading dimension (the LDA parameter 1643 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 1644 1645 Level: intermediate 1646 1647 .keywords: dense, matrix, LAPACK, BLAS 1648 1649 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation() 1650 @*/ 1651 PetscErrorCode MatSeqDenseSetLDA(Mat B,int lda) 1652 { 1653 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1654 PetscFunctionBegin; 1655 if (lda < B->m) SETERRQ2(PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->m); 1656 b->lda = lda; 1657 PetscFunctionReturn(0); 1658 } 1659 1660 /*MC 1661 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 1662 1663 Options Database Keys: 1664 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 1665 1666 Level: beginner 1667 1668 .seealso: MatCreateSeqDense 1669 M*/ 1670 1671 EXTERN_C_BEGIN 1672 #undef __FUNCT__ 1673 #define __FUNCT__ "MatCreate_SeqDense" 1674 PetscErrorCode MatCreate_SeqDense(Mat B) 1675 { 1676 Mat_SeqDense *b; 1677 PetscErrorCode ierr; 1678 PetscMPIInt size; 1679 1680 PetscFunctionBegin; 1681 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1682 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1683 1684 B->m = B->M = PetscMax(B->m,B->M); 1685 B->n = B->N = PetscMax(B->n,B->N); 1686 1687 ierr = PetscNew(Mat_SeqDense,&b);CHKERRQ(ierr); 1688 ierr = PetscMemzero(b,sizeof(Mat_SeqDense));CHKERRQ(ierr); 1689 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1690 B->factor = 0; 1691 B->mapping = 0; 1692 PetscLogObjectMemory(B,sizeof(struct _p_Mat)); 1693 B->data = (void*)b; 1694 1695 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 1696 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 1697 1698 b->pivots = 0; 1699 b->roworiented = PETSC_TRUE; 1700 b->v = 0; 1701 b->lda = B->m; 1702 1703 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 1704 "MatSeqDenseSetPreallocation_SeqDense", 1705 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 1706 PetscFunctionReturn(0); 1707 } 1708 EXTERN_C_END 1709