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