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