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