1 /*$Id: dense.c,v 1.196 2001/03/23 22:04:44 bsmith Exp balay $*/ 2 /* 3 Defines the basic matrix operations for sequential dense. 4 */ 5 6 #include "src/mat/impls/dense/seq/dense.h" 7 #include "petscblaslapack.h" 8 9 #undef __FUNCT__ 10 #define __FUNCT__ "MatAXPY_SeqDense" 11 int MatAXPY_SeqDense(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 __FUNCT__ 23 #define __FUNCT__ "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 __FUNCT__ 52 #define __FUNCT__ "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 __FUNCT__ 69 #define __FUNCT__ "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 __FUNCT__ 90 #define __FUNCT__ "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 __FUNCT__ 109 #define __FUNCT__ "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 __FUNCT__ 120 #define __FUNCT__ "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 __FUNCT__ 135 #define __FUNCT__ "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 __FUNCT__ 146 #define __FUNCT__ "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 __FUNCT__ 157 #define __FUNCT__ "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 __FUNCT__ 178 #define __FUNCT__ "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 __FUNCT__ 204 #define __FUNCT__ "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 __FUNCT__ 230 #define __FUNCT__ "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 __FUNCT__ 264 #define __FUNCT__ "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 __FUNCT__ 302 #define __FUNCT__ "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 __FUNCT__ 363 #define __FUNCT__ "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 __FUNCT__ 382 #define __FUNCT__ "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 __FUNCT__ 401 #define __FUNCT__ "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 __FUNCT__ 421 #define __FUNCT__ "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 __FUNCT__ 443 #define __FUNCT__ "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 __FUNCT__ 465 #define __FUNCT__ "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 __FUNCT__ 476 #define __FUNCT__ "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 __FUNCT__ 549 #define __FUNCT__ "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 __FUNCT__ 572 #define __FUNCT__ "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 __FUNCT__ 647 #define __FUNCT__ "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)A,&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 __FUNCT__ 725 #define __FUNCT__ "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 __FUNCT__ 793 #define __FUNCT__ "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 __FUNCT__ 865 #define __FUNCT__ "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 __FUNCT__ 888 #define __FUNCT__ "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 __FUNCT__ 916 #define __FUNCT__ "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 __FUNCT__ 933 #define __FUNCT__ "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 __FUNCT__ 980 #define __FUNCT__ "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 __FUNCT__ 1003 #define __FUNCT__ "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 __FUNCT__ 1024 #define __FUNCT__ "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 __FUNCT__ 1060 #define __FUNCT__ "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 __FUNCT__ 1107 #define __FUNCT__ "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_NO_NEW_NONZERO_LOCATIONS || 1120 op == MAT_YES_NEW_NONZERO_LOCATIONS || 1121 op == MAT_NEW_NONZERO_LOCATION_ERR || 1122 op == MAT_NO_NEW_DIAGONALS || 1123 op == MAT_YES_NEW_DIAGONALS || 1124 op == MAT_IGNORE_OFF_PROC_ENTRIES || 1125 op == MAT_USE_HASH_TABLE) { 1126 PetscLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1127 } else { 1128 SETERRQ(PETSC_ERR_SUP,"unknown option"); 1129 } 1130 PetscFunctionReturn(0); 1131 } 1132 1133 #undef __FUNCT__ 1134 #define __FUNCT__ "MatZeroEntries_SeqDense" 1135 int MatZeroEntries_SeqDense(Mat A) 1136 { 1137 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1138 int ierr; 1139 1140 PetscFunctionBegin; 1141 ierr = PetscMemzero(l->v,A->m*A->n*sizeof(Scalar));CHKERRQ(ierr); 1142 PetscFunctionReturn(0); 1143 } 1144 1145 #undef __FUNCT__ 1146 #define __FUNCT__ "MatGetBlockSize_SeqDense" 1147 int MatGetBlockSize_SeqDense(Mat A,int *bs) 1148 { 1149 PetscFunctionBegin; 1150 *bs = 1; 1151 PetscFunctionReturn(0); 1152 } 1153 1154 #undef __FUNCT__ 1155 #define __FUNCT__ "MatZeroRows_SeqDense" 1156 int MatZeroRows_SeqDense(Mat A,IS is,Scalar *diag) 1157 { 1158 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1159 int n = A->n,i,j,ierr,N,*rows; 1160 Scalar *slot; 1161 1162 PetscFunctionBegin; 1163 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1164 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1165 for (i=0; i<N; i++) { 1166 slot = l->v + rows[i]; 1167 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1168 } 1169 if (diag) { 1170 for (i=0; i<N; i++) { 1171 slot = l->v + (n+1)*rows[i]; 1172 *slot = *diag; 1173 } 1174 } 1175 ISRestoreIndices(is,&rows); 1176 PetscFunctionReturn(0); 1177 } 1178 1179 #undef __FUNCT__ 1180 #define __FUNCT__ "MatGetOwnershipRange_SeqDense" 1181 int MatGetOwnershipRange_SeqDense(Mat A,int *m,int *n) 1182 { 1183 PetscFunctionBegin; 1184 if (m) *m = 0; 1185 if (n) *n = A->m; 1186 PetscFunctionReturn(0); 1187 } 1188 1189 #undef __FUNCT__ 1190 #define __FUNCT__ "MatGetArray_SeqDense" 1191 int MatGetArray_SeqDense(Mat A,Scalar **array) 1192 { 1193 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1194 1195 PetscFunctionBegin; 1196 *array = mat->v; 1197 PetscFunctionReturn(0); 1198 } 1199 1200 #undef __FUNCT__ 1201 #define __FUNCT__ "MatRestoreArray_SeqDense" 1202 int MatRestoreArray_SeqDense(Mat A,Scalar **array) 1203 { 1204 PetscFunctionBegin; 1205 *array = 0; /* user cannot accidently use the array later */ 1206 PetscFunctionReturn(0); 1207 } 1208 1209 #undef __FUNCT__ 1210 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1211 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatReuse scall,Mat *B) 1212 { 1213 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1214 int i,j,ierr,m = A->m,*irow,*icol,nrows,ncols; 1215 Scalar *av,*bv,*v = mat->v; 1216 Mat newmat; 1217 1218 PetscFunctionBegin; 1219 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1220 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1221 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1222 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1223 1224 /* Check submatrixcall */ 1225 if (scall == MAT_REUSE_MATRIX) { 1226 int n_cols,n_rows; 1227 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1228 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1229 newmat = *B; 1230 } else { 1231 /* Create and fill new matrix */ 1232 ierr = MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat);CHKERRQ(ierr); 1233 } 1234 1235 /* Now extract the data pointers and do the copy,column at a time */ 1236 bv = ((Mat_SeqDense*)newmat->data)->v; 1237 1238 for (i=0; i<ncols; i++) { 1239 av = v + m*icol[i]; 1240 for (j=0; j<nrows; j++) { 1241 *bv++ = av[irow[j]]; 1242 } 1243 } 1244 1245 /* Assemble the matrices so that the correct flags are set */ 1246 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1247 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1248 1249 /* Free work space */ 1250 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1251 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1252 *B = newmat; 1253 PetscFunctionReturn(0); 1254 } 1255 1256 #undef __FUNCT__ 1257 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1258 int MatGetSubMatrices_SeqDense(Mat A,int n,IS *irow,IS *icol,MatReuse scall,Mat **B) 1259 { 1260 int ierr,i; 1261 1262 PetscFunctionBegin; 1263 if (scall == MAT_INITIAL_MATRIX) { 1264 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1265 } 1266 1267 for (i=0; i<n; i++) { 1268 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1269 } 1270 PetscFunctionReturn(0); 1271 } 1272 1273 #undef __FUNCT__ 1274 #define __FUNCT__ "MatCopy_SeqDense" 1275 int MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1276 { 1277 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1278 int ierr; 1279 PetscTruth flag; 1280 1281 PetscFunctionBegin; 1282 ierr = PetscTypeCompare((PetscObject)B,MATSEQDENSE,&flag);CHKERRQ(ierr); 1283 if (!flag) { 1284 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1285 PetscFunctionReturn(0); 1286 } 1287 if (A->m != B->m || A->n != B->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1288 ierr = PetscMemcpy(b->v,a->v,A->m*A->n*sizeof(Scalar));CHKERRQ(ierr); 1289 PetscFunctionReturn(0); 1290 } 1291 1292 #undef __FUNCT__ 1293 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1294 int MatSetUpPreallocation_SeqDense(Mat A) 1295 { 1296 int ierr; 1297 1298 PetscFunctionBegin; 1299 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1300 PetscFunctionReturn(0); 1301 } 1302 1303 /* -------------------------------------------------------------------*/ 1304 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1305 MatGetRow_SeqDense, 1306 MatRestoreRow_SeqDense, 1307 MatMult_SeqDense, 1308 MatMultAdd_SeqDense, 1309 MatMultTranspose_SeqDense, 1310 MatMultTransposeAdd_SeqDense, 1311 MatSolve_SeqDense, 1312 MatSolveAdd_SeqDense, 1313 MatSolveTranspose_SeqDense, 1314 MatSolveTransposeAdd_SeqDense, 1315 MatLUFactor_SeqDense, 1316 MatCholeskyFactor_SeqDense, 1317 MatRelax_SeqDense, 1318 MatTranspose_SeqDense, 1319 MatGetInfo_SeqDense, 1320 MatEqual_SeqDense, 1321 MatGetDiagonal_SeqDense, 1322 MatDiagonalScale_SeqDense, 1323 MatNorm_SeqDense, 1324 0, 1325 0, 1326 0, 1327 MatSetOption_SeqDense, 1328 MatZeroEntries_SeqDense, 1329 MatZeroRows_SeqDense, 1330 MatLUFactorSymbolic_SeqDense, 1331 MatLUFactorNumeric_SeqDense, 1332 MatCholeskyFactorSymbolic_SeqDense, 1333 MatCholeskyFactorNumeric_SeqDense, 1334 MatSetUpPreallocation_SeqDense, 1335 0, 1336 MatGetOwnershipRange_SeqDense, 1337 0, 1338 0, 1339 MatGetArray_SeqDense, 1340 MatRestoreArray_SeqDense, 1341 MatDuplicate_SeqDense, 1342 0, 1343 0, 1344 0, 1345 0, 1346 MatAXPY_SeqDense, 1347 MatGetSubMatrices_SeqDense, 1348 0, 1349 MatGetValues_SeqDense, 1350 MatCopy_SeqDense, 1351 0, 1352 MatScale_SeqDense, 1353 0, 1354 0, 1355 0, 1356 MatGetBlockSize_SeqDense, 1357 0, 1358 0, 1359 0, 1360 0, 1361 0, 1362 0, 1363 0, 1364 0, 1365 0, 1366 0, 1367 MatDestroy_SeqDense, 1368 MatView_SeqDense, 1369 MatGetMaps_Petsc}; 1370 1371 #undef __FUNCT__ 1372 #define __FUNCT__ "MatCreateSeqDense" 1373 /*@C 1374 MatCreateSeqDense - Creates a sequential dense matrix that 1375 is stored in column major order (the usual Fortran 77 manner). Many 1376 of the matrix operations use the BLAS and LAPACK routines. 1377 1378 Collective on MPI_Comm 1379 1380 Input Parameters: 1381 + comm - MPI communicator, set to PETSC_COMM_SELF 1382 . m - number of rows 1383 . n - number of columns 1384 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1385 to control all matrix memory allocation. 1386 1387 Output Parameter: 1388 . A - the matrix 1389 1390 Notes: 1391 The data input variable is intended primarily for Fortran programmers 1392 who wish to allocate their own matrix memory space. Most users should 1393 set data=PETSC_NULL. 1394 1395 Level: intermediate 1396 1397 .keywords: dense, matrix, LAPACK, BLAS 1398 1399 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1400 @*/ 1401 int MatCreateSeqDense(MPI_Comm comm,int m,int n,Scalar *data,Mat *A) 1402 { 1403 int ierr; 1404 1405 PetscFunctionBegin; 1406 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1407 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1408 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1409 PetscFunctionReturn(0); 1410 } 1411 1412 #undef __FUNCT__ 1413 #define __FUNCT__ "MatSeqDensePreallocation" 1414 /*@C 1415 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1416 1417 Collective on MPI_Comm 1418 1419 Input Parameters: 1420 + A - the matrix 1421 - data - the array (or PETSC_NULL) 1422 1423 Notes: 1424 The data input variable is intended primarily for Fortran programmers 1425 who wish to allocate their own matrix memory space. Most users should 1426 set data=PETSC_NULL. 1427 1428 Level: intermediate 1429 1430 .keywords: dense, matrix, LAPACK, BLAS 1431 1432 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1433 @*/ 1434 int MatSeqDenseSetPreallocation(Mat B,Scalar *data) 1435 { 1436 Mat_SeqDense *b; 1437 int ierr; 1438 PetscTruth flg2; 1439 1440 PetscFunctionBegin; 1441 ierr = PetscTypeCompare((PetscObject)B,MATSEQDENSE,&flg2);CHKERRQ(ierr); 1442 if (!flg2) PetscFunctionReturn(0); 1443 B->preallocated = PETSC_TRUE; 1444 b = (Mat_SeqDense*)B->data; 1445 if (!data) { 1446 ierr = PetscMalloc((B->m*B->n+1)*sizeof(Scalar),&b->v);CHKERRQ(ierr); 1447 ierr = PetscMemzero(b->v,B->m*B->n*sizeof(Scalar));CHKERRQ(ierr); 1448 b->user_alloc = PETSC_FALSE; 1449 PetscLogObjectMemory(B,B->n*B->m*sizeof(Scalar)); 1450 } else { /* user-allocated storage */ 1451 b->v = data; 1452 b->user_alloc = PETSC_TRUE; 1453 } 1454 PetscFunctionReturn(0); 1455 } 1456 1457 EXTERN_C_BEGIN 1458 #undef __FUNCT__ 1459 #define __FUNCT__ "MatCreate_SeqDense" 1460 int MatCreate_SeqDense(Mat B) 1461 { 1462 Mat_SeqDense *b; 1463 int ierr,size; 1464 1465 PetscFunctionBegin; 1466 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1467 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1468 1469 B->m = B->M = PetscMax(B->m,B->M); 1470 B->n = B->N = PetscMax(B->n,B->N); 1471 1472 ierr = PetscNew(Mat_SeqDense,&b);CHKERRQ(ierr); 1473 ierr = PetscMemzero(b,sizeof(Mat_SeqDense));CHKERRQ(ierr); 1474 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1475 B->factor = 0; 1476 B->mapping = 0; 1477 PetscLogObjectMemory(B,sizeof(struct _p_Mat)); 1478 B->data = (void*)b; 1479 1480 ierr = MapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 1481 ierr = MapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 1482 1483 b->pivots = 0; 1484 b->roworiented = PETSC_TRUE; 1485 b->v = 0; 1486 1487 PetscFunctionReturn(0); 1488 } 1489 EXTERN_C_END 1490 1491 1492 1493 1494