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