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; i++) { 765 v = a->v + i; 766 for (j=0; j<A->n; j++) { 767 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;} 768 v += a->lda; 769 } 770 } 771 #endif 772 if (format == PETSC_VIEWER_ASCII_MATLAB) { 773 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 774 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %d %d \n",A->m,A->n);CHKERRQ(ierr); 775 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%d,%d);\n",name,A->m,A->n);CHKERRQ(ierr); 776 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 777 } 778 779 for (i=0; i<A->m; i++) { 780 v = a->v + i; 781 for (j=0; j<A->n; j++) { 782 #if defined(PETSC_USE_COMPLEX) 783 if (allreal) { 784 ierr = PetscViewerASCIIPrintf(viewer,"%6.4e ",PetscRealPart(*v));CHKERRQ(ierr); 785 } else { 786 ierr = PetscViewerASCIIPrintf(viewer,"%6.4e + %6.4e i ",PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 787 } 788 #else 789 ierr = PetscViewerASCIIPrintf(viewer,"%6.4e ",*v);CHKERRQ(ierr); 790 #endif 791 v += a->lda; 792 } 793 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 794 } 795 if (format == PETSC_VIEWER_ASCII_MATLAB) { 796 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 797 } 798 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 799 } 800 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 801 PetscFunctionReturn(0); 802 } 803 804 #undef __FUNCT__ 805 #define __FUNCT__ "MatView_SeqDense_Binary" 806 static int MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 807 { 808 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 809 int ict,j,n = A->n,m = A->m,i,fd,*col_lens,ierr,nz = m*n; 810 PetscScalar *v,*anonz,*vals; 811 PetscViewerFormat format; 812 813 PetscFunctionBegin; 814 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 815 816 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 817 if (format == PETSC_VIEWER_BINARY_NATIVE) { 818 /* store the matrix as a dense matrix */ 819 ierr = PetscMalloc(4*sizeof(int),&col_lens);CHKERRQ(ierr); 820 col_lens[0] = MAT_FILE_COOKIE; 821 col_lens[1] = m; 822 col_lens[2] = n; 823 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 824 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,1);CHKERRQ(ierr); 825 ierr = PetscFree(col_lens);CHKERRQ(ierr); 826 827 /* write out matrix, by rows */ 828 ierr = PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 829 v = a->v; 830 for (i=0; i<m; i++) { 831 for (j=0; j<n; j++) { 832 vals[i + j*m] = *v++; 833 } 834 } 835 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,0);CHKERRQ(ierr); 836 ierr = PetscFree(vals);CHKERRQ(ierr); 837 } else { 838 ierr = PetscMalloc((4+nz)*sizeof(int),&col_lens);CHKERRQ(ierr); 839 col_lens[0] = MAT_FILE_COOKIE; 840 col_lens[1] = m; 841 col_lens[2] = n; 842 col_lens[3] = nz; 843 844 /* store lengths of each row and write (including header) to file */ 845 for (i=0; i<m; i++) col_lens[4+i] = n; 846 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,1);CHKERRQ(ierr); 847 848 /* Possibly should write in smaller increments, not whole matrix at once? */ 849 /* store column indices (zero start index) */ 850 ict = 0; 851 for (i=0; i<m; i++) { 852 for (j=0; j<n; j++) col_lens[ict++] = j; 853 } 854 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,0);CHKERRQ(ierr); 855 ierr = PetscFree(col_lens);CHKERRQ(ierr); 856 857 /* store nonzero values */ 858 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);CHKERRQ(ierr); 859 ict = 0; 860 for (i=0; i<m; i++) { 861 v = a->v + i; 862 for (j=0; j<n; j++) { 863 anonz[ict++] = *v; v += a->lda; 864 } 865 } 866 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,0);CHKERRQ(ierr); 867 ierr = PetscFree(anonz);CHKERRQ(ierr); 868 } 869 PetscFunctionReturn(0); 870 } 871 872 #undef __FUNCT__ 873 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 874 int MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 875 { 876 Mat A = (Mat) Aa; 877 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 878 int m = A->m,n = A->n,color,i,j,ierr; 879 PetscScalar *v = a->v; 880 PetscViewer viewer; 881 PetscDraw popup; 882 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 883 PetscViewerFormat format; 884 885 PetscFunctionBegin; 886 887 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 888 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 889 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 890 891 /* Loop over matrix elements drawing boxes */ 892 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 893 /* Blue for negative and Red for positive */ 894 color = PETSC_DRAW_BLUE; 895 for(j = 0; j < n; j++) { 896 x_l = j; 897 x_r = x_l + 1.0; 898 for(i = 0; i < m; i++) { 899 y_l = m - i - 1.0; 900 y_r = y_l + 1.0; 901 #if defined(PETSC_USE_COMPLEX) 902 if (PetscRealPart(v[j*m+i]) > 0.) { 903 color = PETSC_DRAW_RED; 904 } else if (PetscRealPart(v[j*m+i]) < 0.) { 905 color = PETSC_DRAW_BLUE; 906 } else { 907 continue; 908 } 909 #else 910 if (v[j*m+i] > 0.) { 911 color = PETSC_DRAW_RED; 912 } else if (v[j*m+i] < 0.) { 913 color = PETSC_DRAW_BLUE; 914 } else { 915 continue; 916 } 917 #endif 918 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 919 } 920 } 921 } else { 922 /* use contour shading to indicate magnitude of values */ 923 /* first determine max of all nonzero values */ 924 for(i = 0; i < m*n; i++) { 925 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 926 } 927 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 928 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 929 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 930 for(j = 0; j < n; j++) { 931 x_l = j; 932 x_r = x_l + 1.0; 933 for(i = 0; i < m; i++) { 934 y_l = m - i - 1.0; 935 y_r = y_l + 1.0; 936 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 937 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 938 } 939 } 940 } 941 PetscFunctionReturn(0); 942 } 943 944 #undef __FUNCT__ 945 #define __FUNCT__ "MatView_SeqDense_Draw" 946 int MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 947 { 948 PetscDraw draw; 949 PetscTruth isnull; 950 PetscReal xr,yr,xl,yl,h,w; 951 int ierr; 952 953 PetscFunctionBegin; 954 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 955 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 956 if (isnull == PETSC_TRUE) PetscFunctionReturn(0); 957 958 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 959 xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0; 960 xr += w; yr += h; xl = -w; yl = -h; 961 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 962 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 963 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 964 PetscFunctionReturn(0); 965 } 966 967 #undef __FUNCT__ 968 #define __FUNCT__ "MatView_SeqDense" 969 int MatView_SeqDense(Mat A,PetscViewer viewer) 970 { 971 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 972 int ierr; 973 PetscTruth issocket,isascii,isbinary,isdraw; 974 975 PetscFunctionBegin; 976 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 977 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 978 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 979 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 980 981 if (issocket) { 982 if (a->lda>A->m) SETERRQ(1,"Case can not handle LDA"); 983 ierr = PetscViewerSocketPutScalar(viewer,A->m,A->n,a->v);CHKERRQ(ierr); 984 } else if (isascii) { 985 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 986 } else if (isbinary) { 987 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 988 } else if (isdraw) { 989 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 990 } else { 991 SETERRQ1(1,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name); 992 } 993 PetscFunctionReturn(0); 994 } 995 996 #undef __FUNCT__ 997 #define __FUNCT__ "MatDestroy_SeqDense" 998 int MatDestroy_SeqDense(Mat mat) 999 { 1000 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1001 int ierr; 1002 1003 PetscFunctionBegin; 1004 #if defined(PETSC_USE_LOG) 1005 PetscLogObjectState((PetscObject)mat,"Rows %d Cols %d",mat->m,mat->n); 1006 #endif 1007 if (l->pivots) {ierr = PetscFree(l->pivots);CHKERRQ(ierr);} 1008 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1009 ierr = PetscFree(l);CHKERRQ(ierr); 1010 PetscFunctionReturn(0); 1011 } 1012 1013 #undef __FUNCT__ 1014 #define __FUNCT__ "MatTranspose_SeqDense" 1015 int MatTranspose_SeqDense(Mat A,Mat *matout) 1016 { 1017 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1018 int k,j,m,n,M,ierr; 1019 PetscScalar *v,tmp; 1020 1021 PetscFunctionBegin; 1022 v = mat->v; m = A->m; M = mat->lda; n = A->n; 1023 if (!matout) { /* in place transpose */ 1024 if (m != n) { 1025 SETERRQ(1,"Can not transpose non-square matrix in place"); 1026 } else { 1027 for (j=0; j<m; j++) { 1028 for (k=0; k<j; k++) { 1029 tmp = v[j + k*M]; 1030 v[j + k*M] = v[k + j*M]; 1031 v[k + j*M] = tmp; 1032 } 1033 } 1034 } 1035 } else { /* out-of-place transpose */ 1036 Mat tmat; 1037 Mat_SeqDense *tmatd; 1038 PetscScalar *v2; 1039 1040 ierr = MatCreateSeqDense(A->comm,A->n,A->m,PETSC_NULL,&tmat);CHKERRQ(ierr); 1041 tmatd = (Mat_SeqDense*)tmat->data; 1042 v = mat->v; v2 = tmatd->v; 1043 for (j=0; j<n; j++) { 1044 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1045 } 1046 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1047 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1048 *matout = tmat; 1049 } 1050 PetscFunctionReturn(0); 1051 } 1052 1053 #undef __FUNCT__ 1054 #define __FUNCT__ "MatEqual_SeqDense" 1055 int MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg) 1056 { 1057 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1058 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1059 int i,j; 1060 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1061 1062 PetscFunctionBegin; 1063 if (A1->m != A2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1064 if (A1->n != A2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1065 for (i=0; i<A1->m; i++) { 1066 v1 = mat1->v+i; v2 = mat2->v+i; 1067 for (j=0; j<A1->n; j++) { 1068 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1069 v1 += mat1->lda; v2 += mat2->lda; 1070 } 1071 } 1072 *flg = PETSC_TRUE; 1073 PetscFunctionReturn(0); 1074 } 1075 1076 #undef __FUNCT__ 1077 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1078 int MatGetDiagonal_SeqDense(Mat A,Vec v) 1079 { 1080 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1081 int ierr,i,n,len; 1082 PetscScalar *x,zero = 0.0; 1083 1084 PetscFunctionBegin; 1085 ierr = VecSet(&zero,v);CHKERRQ(ierr); 1086 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1087 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1088 len = PetscMin(A->m,A->n); 1089 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1090 for (i=0; i<len; i++) { 1091 x[i] = mat->v[i*mat->lda + i]; 1092 } 1093 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1094 PetscFunctionReturn(0); 1095 } 1096 1097 #undef __FUNCT__ 1098 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1099 int MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1100 { 1101 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1102 PetscScalar *l,*r,x,*v; 1103 int ierr,i,j,m = A->m,n = A->n; 1104 1105 PetscFunctionBegin; 1106 if (ll) { 1107 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1108 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1109 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1110 for (i=0; i<m; i++) { 1111 x = l[i]; 1112 v = mat->v + i; 1113 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1114 } 1115 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1116 PetscLogFlops(n*m); 1117 } 1118 if (rr) { 1119 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1120 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1121 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1122 for (i=0; i<n; i++) { 1123 x = r[i]; 1124 v = mat->v + i*m; 1125 for (j=0; j<m; j++) { (*v++) *= x;} 1126 } 1127 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1128 PetscLogFlops(n*m); 1129 } 1130 PetscFunctionReturn(0); 1131 } 1132 1133 #undef __FUNCT__ 1134 #define __FUNCT__ "MatNorm_SeqDense" 1135 int MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1136 { 1137 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1138 PetscScalar *v = mat->v; 1139 PetscReal sum = 0.0; 1140 int lda=mat->lda,m=A->m,i,j; 1141 1142 PetscFunctionBegin; 1143 if (type == NORM_FROBENIUS) { 1144 if (lda>m) { 1145 for (j=0; j<A->n; j++) { 1146 v = mat->v+j*lda; 1147 for (i=0; i<m; i++) { 1148 #if defined(PETSC_USE_COMPLEX) 1149 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1150 #else 1151 sum += (*v)*(*v); v++; 1152 #endif 1153 } 1154 } 1155 } else { 1156 for (i=0; i<A->n*A->m; i++) { 1157 #if defined(PETSC_USE_COMPLEX) 1158 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1159 #else 1160 sum += (*v)*(*v); v++; 1161 #endif 1162 } 1163 } 1164 *nrm = sqrt(sum); 1165 PetscLogFlops(2*A->n*A->m); 1166 } else if (type == NORM_1) { 1167 *nrm = 0.0; 1168 for (j=0; j<A->n; j++) { 1169 v = mat->v + j*mat->lda; 1170 sum = 0.0; 1171 for (i=0; i<A->m; i++) { 1172 sum += PetscAbsScalar(*v); v++; 1173 } 1174 if (sum > *nrm) *nrm = sum; 1175 } 1176 PetscLogFlops(A->n*A->m); 1177 } else if (type == NORM_INFINITY) { 1178 *nrm = 0.0; 1179 for (j=0; j<A->m; j++) { 1180 v = mat->v + j; 1181 sum = 0.0; 1182 for (i=0; i<A->n; i++) { 1183 sum += PetscAbsScalar(*v); v += mat->lda; 1184 } 1185 if (sum > *nrm) *nrm = sum; 1186 } 1187 PetscLogFlops(A->n*A->m); 1188 } else { 1189 SETERRQ(PETSC_ERR_SUP,"No two norm"); 1190 } 1191 PetscFunctionReturn(0); 1192 } 1193 1194 #undef __FUNCT__ 1195 #define __FUNCT__ "MatSetOption_SeqDense" 1196 int MatSetOption_SeqDense(Mat A,MatOption op) 1197 { 1198 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1199 1200 PetscFunctionBegin; 1201 switch (op) { 1202 case MAT_ROW_ORIENTED: 1203 aij->roworiented = PETSC_TRUE; 1204 break; 1205 case MAT_COLUMN_ORIENTED: 1206 aij->roworiented = PETSC_FALSE; 1207 break; 1208 case MAT_ROWS_SORTED: 1209 case MAT_ROWS_UNSORTED: 1210 case MAT_COLUMNS_SORTED: 1211 case MAT_COLUMNS_UNSORTED: 1212 case MAT_NO_NEW_NONZERO_LOCATIONS: 1213 case MAT_YES_NEW_NONZERO_LOCATIONS: 1214 case MAT_NEW_NONZERO_LOCATION_ERR: 1215 case MAT_NO_NEW_DIAGONALS: 1216 case MAT_YES_NEW_DIAGONALS: 1217 case MAT_IGNORE_OFF_PROC_ENTRIES: 1218 case MAT_USE_HASH_TABLE: 1219 PetscLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1220 break; 1221 default: 1222 SETERRQ(PETSC_ERR_SUP,"unknown option"); 1223 } 1224 PetscFunctionReturn(0); 1225 } 1226 1227 #undef __FUNCT__ 1228 #define __FUNCT__ "MatZeroEntries_SeqDense" 1229 int MatZeroEntries_SeqDense(Mat A) 1230 { 1231 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1232 int lda=l->lda,m=A->m,j, ierr; 1233 1234 PetscFunctionBegin; 1235 if (lda>m) { 1236 for (j=0; j<A->n; j++) { 1237 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar)); CHKERRQ(ierr); 1238 } 1239 } else { 1240 ierr = PetscMemzero(l->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 1241 } 1242 PetscFunctionReturn(0); 1243 } 1244 1245 #undef __FUNCT__ 1246 #define __FUNCT__ "MatGetBlockSize_SeqDense" 1247 int MatGetBlockSize_SeqDense(Mat A,int *bs) 1248 { 1249 PetscFunctionBegin; 1250 *bs = 1; 1251 PetscFunctionReturn(0); 1252 } 1253 1254 #undef __FUNCT__ 1255 #define __FUNCT__ "MatZeroRows_SeqDense" 1256 int MatZeroRows_SeqDense(Mat A,IS is,const PetscScalar *diag) 1257 { 1258 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1259 int n = A->n,i,j,ierr,N,*rows; 1260 PetscScalar *slot; 1261 1262 PetscFunctionBegin; 1263 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1264 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1265 for (i=0; i<N; i++) { 1266 slot = l->v + rows[i]; 1267 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1268 } 1269 if (diag) { 1270 for (i=0; i<N; i++) { 1271 slot = l->v + (n+1)*rows[i]; 1272 *slot = *diag; 1273 } 1274 } 1275 ISRestoreIndices(is,&rows); 1276 PetscFunctionReturn(0); 1277 } 1278 1279 #undef __FUNCT__ 1280 #define __FUNCT__ "MatGetArray_SeqDense" 1281 int MatGetArray_SeqDense(Mat A,PetscScalar *array[]) 1282 { 1283 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1284 1285 PetscFunctionBegin; 1286 *array = mat->v; 1287 PetscFunctionReturn(0); 1288 } 1289 1290 #undef __FUNCT__ 1291 #define __FUNCT__ "MatRestoreArray_SeqDense" 1292 int MatRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1293 { 1294 PetscFunctionBegin; 1295 *array = 0; /* user cannot accidently use the array later */ 1296 PetscFunctionReturn(0); 1297 } 1298 1299 #undef __FUNCT__ 1300 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1301 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatReuse scall,Mat *B) 1302 { 1303 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1304 int i,j,ierr,m = A->m,*irow,*icol,nrows,ncols; 1305 PetscScalar *av,*bv,*v = mat->v; 1306 Mat newmat; 1307 1308 PetscFunctionBegin; 1309 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1310 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1311 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1312 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1313 1314 /* Check submatrixcall */ 1315 if (scall == MAT_REUSE_MATRIX) { 1316 int n_cols,n_rows; 1317 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1318 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1319 newmat = *B; 1320 } else { 1321 /* Create and fill new matrix */ 1322 ierr = MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat);CHKERRQ(ierr); 1323 } 1324 1325 /* Now extract the data pointers and do the copy,column at a time */ 1326 bv = ((Mat_SeqDense*)newmat->data)->v; 1327 1328 for (i=0; i<ncols; i++) { 1329 av = v + m*icol[i]; 1330 for (j=0; j<nrows; j++) { 1331 *bv++ = av[irow[j]]; 1332 } 1333 } 1334 1335 /* Assemble the matrices so that the correct flags are set */ 1336 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1337 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1338 1339 /* Free work space */ 1340 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1341 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1342 *B = newmat; 1343 PetscFunctionReturn(0); 1344 } 1345 1346 #undef __FUNCT__ 1347 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1348 int MatGetSubMatrices_SeqDense(Mat A,int n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1349 { 1350 int ierr,i; 1351 1352 PetscFunctionBegin; 1353 if (scall == MAT_INITIAL_MATRIX) { 1354 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1355 } 1356 1357 for (i=0; i<n; i++) { 1358 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1359 } 1360 PetscFunctionReturn(0); 1361 } 1362 1363 #undef __FUNCT__ 1364 #define __FUNCT__ "MatCopy_SeqDense" 1365 int MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1366 { 1367 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1368 int lda1=a->lda,lda2=b->lda, m=A->m,n=A->n, j,ierr; 1369 1370 PetscFunctionBegin; 1371 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1372 if (A->ops->copy != B->ops->copy) { 1373 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1374 PetscFunctionReturn(0); 1375 } 1376 if (m != B->m || n != B->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1377 if (lda1>m || lda2>m) { 1378 for (j=0; j<n; j++) { 1379 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1380 } 1381 } else { 1382 ierr = PetscMemcpy(b->v,a->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 1383 } 1384 PetscFunctionReturn(0); 1385 } 1386 1387 #undef __FUNCT__ 1388 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1389 int MatSetUpPreallocation_SeqDense(Mat A) 1390 { 1391 int ierr; 1392 1393 PetscFunctionBegin; 1394 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1395 PetscFunctionReturn(0); 1396 } 1397 1398 /* -------------------------------------------------------------------*/ 1399 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1400 MatGetRow_SeqDense, 1401 MatRestoreRow_SeqDense, 1402 MatMult_SeqDense, 1403 MatMultAdd_SeqDense, 1404 MatMultTranspose_SeqDense, 1405 MatMultTransposeAdd_SeqDense, 1406 MatSolve_SeqDense, 1407 MatSolveAdd_SeqDense, 1408 MatSolveTranspose_SeqDense, 1409 MatSolveTransposeAdd_SeqDense, 1410 MatLUFactor_SeqDense, 1411 MatCholeskyFactor_SeqDense, 1412 MatRelax_SeqDense, 1413 MatTranspose_SeqDense, 1414 MatGetInfo_SeqDense, 1415 MatEqual_SeqDense, 1416 MatGetDiagonal_SeqDense, 1417 MatDiagonalScale_SeqDense, 1418 MatNorm_SeqDense, 1419 0, 1420 0, 1421 0, 1422 MatSetOption_SeqDense, 1423 MatZeroEntries_SeqDense, 1424 MatZeroRows_SeqDense, 1425 MatLUFactorSymbolic_SeqDense, 1426 MatLUFactorNumeric_SeqDense, 1427 MatCholeskyFactorSymbolic_SeqDense, 1428 MatCholeskyFactorNumeric_SeqDense, 1429 MatSetUpPreallocation_SeqDense, 1430 0, 1431 0, 1432 MatGetArray_SeqDense, 1433 MatRestoreArray_SeqDense, 1434 MatDuplicate_SeqDense, 1435 0, 1436 0, 1437 0, 1438 0, 1439 MatAXPY_SeqDense, 1440 MatGetSubMatrices_SeqDense, 1441 0, 1442 MatGetValues_SeqDense, 1443 MatCopy_SeqDense, 1444 0, 1445 MatScale_SeqDense, 1446 0, 1447 0, 1448 0, 1449 MatGetBlockSize_SeqDense, 1450 0, 1451 0, 1452 0, 1453 0, 1454 0, 1455 0, 1456 0, 1457 0, 1458 0, 1459 0, 1460 MatDestroy_SeqDense, 1461 MatView_SeqDense, 1462 MatGetPetscMaps_Petsc}; 1463 1464 #undef __FUNCT__ 1465 #define __FUNCT__ "MatCreateSeqDense" 1466 /*@C 1467 MatCreateSeqDense - Creates a sequential dense matrix that 1468 is stored in column major order (the usual Fortran 77 manner). Many 1469 of the matrix operations use the BLAS and LAPACK routines. 1470 1471 Collective on MPI_Comm 1472 1473 Input Parameters: 1474 + comm - MPI communicator, set to PETSC_COMM_SELF 1475 . m - number of rows 1476 . n - number of columns 1477 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1478 to control all matrix memory allocation. 1479 1480 Output Parameter: 1481 . A - the matrix 1482 1483 Notes: 1484 The data input variable is intended primarily for Fortran programmers 1485 who wish to allocate their own matrix memory space. Most users should 1486 set data=PETSC_NULL. 1487 1488 Level: intermediate 1489 1490 .keywords: dense, matrix, LAPACK, BLAS 1491 1492 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1493 @*/ 1494 int MatCreateSeqDense(MPI_Comm comm,int m,int n,PetscScalar *data,Mat *A) 1495 { 1496 int ierr; 1497 1498 PetscFunctionBegin; 1499 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1500 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1501 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1502 PetscFunctionReturn(0); 1503 } 1504 1505 #undef __FUNCT__ 1506 #define __FUNCT__ "MatSeqDensePreallocation" 1507 /*@C 1508 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1509 1510 Collective on MPI_Comm 1511 1512 Input Parameters: 1513 + A - the matrix 1514 - data - the array (or PETSC_NULL) 1515 1516 Notes: 1517 The data input variable is intended primarily for Fortran programmers 1518 who wish to allocate their own matrix memory space. Most users should 1519 set data=PETSC_NULL. 1520 1521 Level: intermediate 1522 1523 .keywords: dense, matrix, LAPACK, BLAS 1524 1525 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1526 @*/ 1527 int MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 1528 { 1529 int ierr,(*f)(Mat,PetscScalar[]); 1530 1531 PetscFunctionBegin; 1532 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1533 if (f) { 1534 ierr = (*f)(B,data);CHKERRQ(ierr); 1535 } 1536 PetscFunctionReturn(0); 1537 } 1538 1539 EXTERN_C_BEGIN 1540 #undef __FUNCT__ 1541 #define __FUNCT__ "MatSeqDensePreallocation_SeqDense" 1542 int MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1543 { 1544 Mat_SeqDense *b; 1545 int ierr; 1546 1547 PetscFunctionBegin; 1548 B->preallocated = PETSC_TRUE; 1549 b = (Mat_SeqDense*)B->data; 1550 if (!data) { 1551 ierr = PetscMalloc((B->m*B->n+1)*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1552 ierr = PetscMemzero(b->v,B->m*B->n*sizeof(PetscScalar));CHKERRQ(ierr); 1553 b->user_alloc = PETSC_FALSE; 1554 PetscLogObjectMemory(B,B->n*B->m*sizeof(PetscScalar)); 1555 } else { /* user-allocated storage */ 1556 b->v = data; 1557 b->user_alloc = PETSC_TRUE; 1558 } 1559 PetscFunctionReturn(0); 1560 } 1561 EXTERN_C_END 1562 1563 #undef __FUNCT__ 1564 #define __FUNCT__ "MatSeqDenseSetLDA" 1565 /*@C 1566 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 1567 1568 Input parameter: 1569 + A - the matrix 1570 - lda - the leading dimension 1571 1572 Notes: 1573 This routine is to be used in conjunction with MatSeqDenseSetPreallocation; 1574 it asserts that the preallocation has a leading dimension (the LDA parameter 1575 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 1576 1577 Level: intermediate 1578 1579 .keywords: dense, matrix, LAPACK, BLAS 1580 1581 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation() 1582 @*/ 1583 int MatSeqDenseSetLDA(Mat B,int lda) 1584 { 1585 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1586 PetscFunctionBegin; 1587 if (lda<B->m) SETERRQ(1,"LDA must be at least matrix i dimension"); 1588 b->lda = lda; 1589 PetscFunctionReturn(0); 1590 } 1591 1592 EXTERN_C_BEGIN 1593 #undef __FUNCT__ 1594 #define __FUNCT__ "MatCreate_SeqDense" 1595 int MatCreate_SeqDense(Mat B) 1596 { 1597 Mat_SeqDense *b; 1598 int ierr,size; 1599 1600 PetscFunctionBegin; 1601 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1602 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1603 1604 B->m = B->M = PetscMax(B->m,B->M); 1605 B->n = B->N = PetscMax(B->n,B->N); 1606 1607 ierr = PetscNew(Mat_SeqDense,&b);CHKERRQ(ierr); 1608 ierr = PetscMemzero(b,sizeof(Mat_SeqDense));CHKERRQ(ierr); 1609 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1610 B->factor = 0; 1611 B->mapping = 0; 1612 PetscLogObjectMemory(B,sizeof(struct _p_Mat)); 1613 B->data = (void*)b; 1614 1615 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 1616 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 1617 1618 b->pivots = 0; 1619 b->roworiented = PETSC_TRUE; 1620 b->v = 0; 1621 b->lda = B->m; 1622 1623 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 1624 "MatSeqDenseSetPreallocation_SeqDense", 1625 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 1626 PetscFunctionReturn(0); 1627 } 1628 EXTERN_C_END 1629