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