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