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