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 i,j; 1061 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1062 1063 PetscFunctionBegin; 1064 if (A1->m != A2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1065 if (A1->n != A2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1066 for (i=0; i<A1->m; i++) { 1067 v1 = mat1->v+i; v2 = mat2->v+i; 1068 for (j=0; j<A1->n; j++) { 1069 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1070 v1 += mat1->lda; v2 += mat2->lda; 1071 } 1072 } 1073 *flg = PETSC_TRUE; 1074 PetscFunctionReturn(0); 1075 } 1076 1077 #undef __FUNCT__ 1078 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1079 int MatGetDiagonal_SeqDense(Mat A,Vec v) 1080 { 1081 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1082 int ierr,i,n,len; 1083 PetscScalar *x,zero = 0.0; 1084 1085 PetscFunctionBegin; 1086 ierr = VecSet(&zero,v);CHKERRQ(ierr); 1087 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1088 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1089 len = PetscMin(A->m,A->n); 1090 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1091 for (i=0; i<len; i++) { 1092 x[i] = mat->v[i*mat->lda + i]; 1093 } 1094 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1095 PetscFunctionReturn(0); 1096 } 1097 1098 #undef __FUNCT__ 1099 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1100 int MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1101 { 1102 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1103 PetscScalar *l,*r,x,*v; 1104 int ierr,i,j,m = A->m,n = A->n; 1105 1106 PetscFunctionBegin; 1107 if (ll) { 1108 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1109 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1110 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1111 for (i=0; i<m; i++) { 1112 x = l[i]; 1113 v = mat->v + i; 1114 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1115 } 1116 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1117 PetscLogFlops(n*m); 1118 } 1119 if (rr) { 1120 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1121 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1122 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1123 for (i=0; i<n; i++) { 1124 x = r[i]; 1125 v = mat->v + i*m; 1126 for (j=0; j<m; j++) { (*v++) *= x;} 1127 } 1128 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1129 PetscLogFlops(n*m); 1130 } 1131 PetscFunctionReturn(0); 1132 } 1133 1134 #undef __FUNCT__ 1135 #define __FUNCT__ "MatNorm_SeqDense" 1136 int MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1137 { 1138 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1139 PetscScalar *v = mat->v; 1140 PetscReal sum = 0.0; 1141 int lda=mat->lda,m=A->m,i,j; 1142 1143 PetscFunctionBegin; 1144 if (type == NORM_FROBENIUS) { 1145 if (lda>m) { 1146 for (j=0; j<A->n; j++) { 1147 v = mat->v+j*lda; 1148 for (i=0; i<m; i++) { 1149 #if defined(PETSC_USE_COMPLEX) 1150 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1151 #else 1152 sum += (*v)*(*v); v++; 1153 #endif 1154 } 1155 } 1156 } else { 1157 for (i=0; i<A->n*A->m; i++) { 1158 #if defined(PETSC_USE_COMPLEX) 1159 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1160 #else 1161 sum += (*v)*(*v); v++; 1162 #endif 1163 } 1164 } 1165 *nrm = sqrt(sum); 1166 PetscLogFlops(2*A->n*A->m); 1167 } else if (type == NORM_1) { 1168 *nrm = 0.0; 1169 for (j=0; j<A->n; j++) { 1170 v = mat->v + j*mat->lda; 1171 sum = 0.0; 1172 for (i=0; i<A->m; i++) { 1173 sum += PetscAbsScalar(*v); v++; 1174 } 1175 if (sum > *nrm) *nrm = sum; 1176 } 1177 PetscLogFlops(A->n*A->m); 1178 } else if (type == NORM_INFINITY) { 1179 *nrm = 0.0; 1180 for (j=0; j<A->m; j++) { 1181 v = mat->v + j; 1182 sum = 0.0; 1183 for (i=0; i<A->n; i++) { 1184 sum += PetscAbsScalar(*v); v += mat->lda; 1185 } 1186 if (sum > *nrm) *nrm = sum; 1187 } 1188 PetscLogFlops(A->n*A->m); 1189 } else { 1190 SETERRQ(PETSC_ERR_SUP,"No two norm"); 1191 } 1192 PetscFunctionReturn(0); 1193 } 1194 1195 #undef __FUNCT__ 1196 #define __FUNCT__ "MatSetOption_SeqDense" 1197 int MatSetOption_SeqDense(Mat A,MatOption op) 1198 { 1199 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1200 1201 PetscFunctionBegin; 1202 switch (op) { 1203 case MAT_ROW_ORIENTED: 1204 aij->roworiented = PETSC_TRUE; 1205 break; 1206 case MAT_COLUMN_ORIENTED: 1207 aij->roworiented = PETSC_FALSE; 1208 break; 1209 case MAT_ROWS_SORTED: 1210 case MAT_ROWS_UNSORTED: 1211 case MAT_COLUMNS_SORTED: 1212 case MAT_COLUMNS_UNSORTED: 1213 case MAT_NO_NEW_NONZERO_LOCATIONS: 1214 case MAT_YES_NEW_NONZERO_LOCATIONS: 1215 case MAT_NEW_NONZERO_LOCATION_ERR: 1216 case MAT_NO_NEW_DIAGONALS: 1217 case MAT_YES_NEW_DIAGONALS: 1218 case MAT_IGNORE_OFF_PROC_ENTRIES: 1219 case MAT_USE_HASH_TABLE: 1220 PetscLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1221 break; 1222 default: 1223 SETERRQ(PETSC_ERR_SUP,"unknown option"); 1224 } 1225 PetscFunctionReturn(0); 1226 } 1227 1228 #undef __FUNCT__ 1229 #define __FUNCT__ "MatZeroEntries_SeqDense" 1230 int MatZeroEntries_SeqDense(Mat A) 1231 { 1232 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1233 int lda=l->lda,m=A->m,j, ierr; 1234 1235 PetscFunctionBegin; 1236 if (lda>m) { 1237 for (j=0; j<A->n; j++) { 1238 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar)); CHKERRQ(ierr); 1239 } 1240 } else { 1241 ierr = PetscMemzero(l->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 1242 } 1243 PetscFunctionReturn(0); 1244 } 1245 1246 #undef __FUNCT__ 1247 #define __FUNCT__ "MatGetBlockSize_SeqDense" 1248 int MatGetBlockSize_SeqDense(Mat A,int *bs) 1249 { 1250 PetscFunctionBegin; 1251 *bs = 1; 1252 PetscFunctionReturn(0); 1253 } 1254 1255 #undef __FUNCT__ 1256 #define __FUNCT__ "MatZeroRows_SeqDense" 1257 int MatZeroRows_SeqDense(Mat A,IS is,PetscScalar *diag) 1258 { 1259 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1260 int n = A->n,i,j,ierr,N,*rows; 1261 PetscScalar *slot; 1262 1263 PetscFunctionBegin; 1264 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1265 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1266 for (i=0; i<N; i++) { 1267 slot = l->v + rows[i]; 1268 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1269 } 1270 if (diag) { 1271 for (i=0; i<N; i++) { 1272 slot = l->v + (n+1)*rows[i]; 1273 *slot = *diag; 1274 } 1275 } 1276 ISRestoreIndices(is,&rows); 1277 PetscFunctionReturn(0); 1278 } 1279 1280 #undef __FUNCT__ 1281 #define __FUNCT__ "MatGetArray_SeqDense" 1282 int MatGetArray_SeqDense(Mat A,PetscScalar **array) 1283 { 1284 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1285 1286 PetscFunctionBegin; 1287 *array = mat->v; 1288 PetscFunctionReturn(0); 1289 } 1290 1291 #undef __FUNCT__ 1292 #define __FUNCT__ "MatRestoreArray_SeqDense" 1293 int MatRestoreArray_SeqDense(Mat A,PetscScalar **array) 1294 { 1295 PetscFunctionBegin; 1296 *array = 0; /* user cannot accidently use the array later */ 1297 PetscFunctionReturn(0); 1298 } 1299 1300 #undef __FUNCT__ 1301 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1302 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatReuse scall,Mat *B) 1303 { 1304 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1305 int i,j,ierr,m = A->m,*irow,*icol,nrows,ncols; 1306 PetscScalar *av,*bv,*v = mat->v; 1307 Mat newmat; 1308 1309 PetscFunctionBegin; 1310 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1311 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1312 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1313 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1314 1315 /* Check submatrixcall */ 1316 if (scall == MAT_REUSE_MATRIX) { 1317 int n_cols,n_rows; 1318 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1319 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1320 newmat = *B; 1321 } else { 1322 /* Create and fill new matrix */ 1323 ierr = MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat);CHKERRQ(ierr); 1324 } 1325 1326 /* Now extract the data pointers and do the copy,column at a time */ 1327 bv = ((Mat_SeqDense*)newmat->data)->v; 1328 1329 for (i=0; i<ncols; i++) { 1330 av = v + m*icol[i]; 1331 for (j=0; j<nrows; j++) { 1332 *bv++ = av[irow[j]]; 1333 } 1334 } 1335 1336 /* Assemble the matrices so that the correct flags are set */ 1337 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1338 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1339 1340 /* Free work space */ 1341 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1342 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1343 *B = newmat; 1344 PetscFunctionReturn(0); 1345 } 1346 1347 #undef __FUNCT__ 1348 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1349 int MatGetSubMatrices_SeqDense(Mat A,int n,IS *irow,IS *icol,MatReuse scall,Mat **B) 1350 { 1351 int ierr,i; 1352 1353 PetscFunctionBegin; 1354 if (scall == MAT_INITIAL_MATRIX) { 1355 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1356 } 1357 1358 for (i=0; i<n; i++) { 1359 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1360 } 1361 PetscFunctionReturn(0); 1362 } 1363 1364 #undef __FUNCT__ 1365 #define __FUNCT__ "MatCopy_SeqDense" 1366 int MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1367 { 1368 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1369 int lda1=a->lda,lda2=b->lda, m=A->m,n=A->n, j,ierr; 1370 1371 PetscFunctionBegin; 1372 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1373 if (A->ops->copy != B->ops->copy) { 1374 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1375 PetscFunctionReturn(0); 1376 } 1377 if (m != B->m || n != B->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1378 if (lda1>m || lda2>m) { 1379 for (j=0; j<n; j++) { 1380 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1381 } 1382 } else { 1383 ierr = PetscMemcpy(b->v,a->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 1384 } 1385 PetscFunctionReturn(0); 1386 } 1387 1388 #undef __FUNCT__ 1389 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1390 int MatSetUpPreallocation_SeqDense(Mat A) 1391 { 1392 int ierr; 1393 1394 PetscFunctionBegin; 1395 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1396 PetscFunctionReturn(0); 1397 } 1398 1399 /* -------------------------------------------------------------------*/ 1400 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1401 MatGetRow_SeqDense, 1402 MatRestoreRow_SeqDense, 1403 MatMult_SeqDense, 1404 MatMultAdd_SeqDense, 1405 MatMultTranspose_SeqDense, 1406 MatMultTransposeAdd_SeqDense, 1407 MatSolve_SeqDense, 1408 MatSolveAdd_SeqDense, 1409 MatSolveTranspose_SeqDense, 1410 MatSolveTransposeAdd_SeqDense, 1411 MatLUFactor_SeqDense, 1412 MatCholeskyFactor_SeqDense, 1413 MatRelax_SeqDense, 1414 MatTranspose_SeqDense, 1415 MatGetInfo_SeqDense, 1416 MatEqual_SeqDense, 1417 MatGetDiagonal_SeqDense, 1418 MatDiagonalScale_SeqDense, 1419 MatNorm_SeqDense, 1420 0, 1421 0, 1422 0, 1423 MatSetOption_SeqDense, 1424 MatZeroEntries_SeqDense, 1425 MatZeroRows_SeqDense, 1426 MatLUFactorSymbolic_SeqDense, 1427 MatLUFactorNumeric_SeqDense, 1428 MatCholeskyFactorSymbolic_SeqDense, 1429 MatCholeskyFactorNumeric_SeqDense, 1430 MatSetUpPreallocation_SeqDense, 1431 0, 1432 0, 1433 MatGetArray_SeqDense, 1434 MatRestoreArray_SeqDense, 1435 MatDuplicate_SeqDense, 1436 0, 1437 0, 1438 0, 1439 0, 1440 MatAXPY_SeqDense, 1441 MatGetSubMatrices_SeqDense, 1442 0, 1443 MatGetValues_SeqDense, 1444 MatCopy_SeqDense, 1445 0, 1446 MatScale_SeqDense, 1447 0, 1448 0, 1449 0, 1450 MatGetBlockSize_SeqDense, 1451 0, 1452 0, 1453 0, 1454 0, 1455 0, 1456 0, 1457 0, 1458 0, 1459 0, 1460 0, 1461 MatDestroy_SeqDense, 1462 MatView_SeqDense, 1463 MatGetPetscMaps_Petsc}; 1464 1465 #undef __FUNCT__ 1466 #define __FUNCT__ "MatCreateSeqDense" 1467 /*@C 1468 MatCreateSeqDense - Creates a sequential dense matrix that 1469 is stored in column major order (the usual Fortran 77 manner). Many 1470 of the matrix operations use the BLAS and LAPACK routines. 1471 1472 Collective on MPI_Comm 1473 1474 Input Parameters: 1475 + comm - MPI communicator, set to PETSC_COMM_SELF 1476 . m - number of rows 1477 . n - number of columns 1478 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1479 to control all matrix memory allocation. 1480 1481 Output Parameter: 1482 . A - the matrix 1483 1484 Notes: 1485 The data input variable is intended primarily for Fortran programmers 1486 who wish to allocate their own matrix memory space. Most users should 1487 set data=PETSC_NULL. 1488 1489 Level: intermediate 1490 1491 .keywords: dense, matrix, LAPACK, BLAS 1492 1493 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1494 @*/ 1495 int MatCreateSeqDense(MPI_Comm comm,int m,int n,PetscScalar *data,Mat *A) 1496 { 1497 int ierr; 1498 1499 PetscFunctionBegin; 1500 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1501 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1502 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1503 PetscFunctionReturn(0); 1504 } 1505 1506 #undef __FUNCT__ 1507 #define __FUNCT__ "MatSeqDensePreallocation" 1508 /*@C 1509 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1510 1511 Collective on MPI_Comm 1512 1513 Input Parameters: 1514 + A - the matrix 1515 - data - the array (or PETSC_NULL) 1516 1517 Notes: 1518 The data input variable is intended primarily for Fortran programmers 1519 who wish to allocate their own matrix memory space. Most users should 1520 set data=PETSC_NULL. 1521 1522 Level: intermediate 1523 1524 .keywords: dense, matrix, LAPACK, BLAS 1525 1526 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1527 @*/ 1528 int MatSeqDenseSetPreallocation(Mat B,PetscScalar *data) 1529 { 1530 int ierr,(*f)(Mat,PetscScalar*); 1531 1532 PetscFunctionBegin; 1533 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1534 if (f) { 1535 ierr = (*f)(B,data);CHKERRQ(ierr); 1536 } 1537 PetscFunctionReturn(0); 1538 } 1539 1540 EXTERN_C_BEGIN 1541 #undef __FUNCT__ 1542 #define __FUNCT__ "MatSeqDensePreallocation_SeqDense" 1543 int MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1544 { 1545 Mat_SeqDense *b; 1546 int ierr; 1547 1548 PetscFunctionBegin; 1549 B->preallocated = PETSC_TRUE; 1550 b = (Mat_SeqDense*)B->data; 1551 if (!data) { 1552 ierr = PetscMalloc((B->m*B->n+1)*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1553 ierr = PetscMemzero(b->v,B->m*B->n*sizeof(PetscScalar));CHKERRQ(ierr); 1554 b->user_alloc = PETSC_FALSE; 1555 PetscLogObjectMemory(B,B->n*B->m*sizeof(PetscScalar)); 1556 } else { /* user-allocated storage */ 1557 b->v = data; 1558 b->user_alloc = PETSC_TRUE; 1559 } 1560 PetscFunctionReturn(0); 1561 } 1562 EXTERN_C_END 1563 1564 #undef __FUNCT__ 1565 #define __FUNCT__ "MatSeqDenseSetLDA" 1566 /*@C 1567 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 1568 1569 Input parameter: 1570 + A - the matrix 1571 - lda - the leading dimension 1572 1573 Notes: 1574 This routine is to be used in conjunction with MatSeqDenseSetPreallocation; 1575 it asserts that the preallocation has a leading dimension (the LDA parameter 1576 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 1577 1578 Level: intermediate 1579 1580 .keywords: dense, matrix, LAPACK, BLAS 1581 1582 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation() 1583 @*/ 1584 int MatSeqDenseSetLDA(Mat B,int lda) 1585 { 1586 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1587 PetscFunctionBegin; 1588 if (lda<B->m) SETERRQ(1,"LDA must be at least matrix i dimension"); 1589 b->lda = lda; 1590 PetscFunctionReturn(0); 1591 } 1592 1593 EXTERN_C_BEGIN 1594 #undef __FUNCT__ 1595 #define __FUNCT__ "MatCreate_SeqDense" 1596 int MatCreate_SeqDense(Mat B) 1597 { 1598 Mat_SeqDense *b; 1599 int ierr,size; 1600 1601 PetscFunctionBegin; 1602 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1603 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1604 1605 B->m = B->M = PetscMax(B->m,B->M); 1606 B->n = B->N = PetscMax(B->n,B->N); 1607 1608 ierr = PetscNew(Mat_SeqDense,&b);CHKERRQ(ierr); 1609 ierr = PetscMemzero(b,sizeof(Mat_SeqDense));CHKERRQ(ierr); 1610 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1611 B->factor = 0; 1612 B->mapping = 0; 1613 PetscLogObjectMemory(B,sizeof(struct _p_Mat)); 1614 B->data = (void*)b; 1615 1616 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 1617 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 1618 1619 b->pivots = 0; 1620 b->roworiented = PETSC_TRUE; 1621 b->v = 0; 1622 b->lda = B->m; 1623 1624 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 1625 "MatSeqDenseSetPreallocation_SeqDense", 1626 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 1627 PetscFunctionReturn(0); 1628 } 1629 EXTERN_C_END 1630