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