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