1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: dense.c,v 1.138 1998/03/12 23:18:15 bsmith Exp balay $"; 3 #endif 4 /* 5 Defines the basic matrix operations for sequential dense. 6 */ 7 8 #include "src/mat/impls/dense/seq/dense.h" 9 #include "pinclude/plapack.h" 10 #include "pinclude/pviewer.h" 11 12 #undef __FUNC__ 13 #define __FUNC__ "MatAXPY_SeqDense" 14 int MatAXPY_SeqDense(Scalar *alpha,Mat X,Mat Y) 15 { 16 Mat_SeqDense *x = (Mat_SeqDense*) X->data,*y = (Mat_SeqDense*) Y->data; 17 int N = x->m*x->n, one = 1; 18 19 PetscFunctionBegin; 20 BLaxpy_( &N, alpha, x->v, &one, y->v, &one ); 21 PLogFlops(2*N-1); 22 PetscFunctionReturn(0); 23 } 24 25 #undef __FUNC__ 26 #define __FUNC__ "MatGetInfo_SeqDense" 27 int MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 28 { 29 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 30 int i,N = a->m*a->n,count = 0; 31 Scalar *v = a->v; 32 33 PetscFunctionBegin; 34 for ( i=0; i<N; i++ ) {if (*v != 0.0) count++; v++;} 35 36 info->rows_global = (double)a->m; 37 info->columns_global = (double)a->n; 38 info->rows_local = (double)a->m; 39 info->columns_local = (double)a->n; 40 info->block_size = 1.0; 41 info->nz_allocated = (double)N; 42 info->nz_used = (double)count; 43 info->nz_unneeded = (double)(N-count); 44 info->assemblies = (double)A->num_ass; 45 info->mallocs = 0; 46 info->memory = A->mem; 47 info->fill_ratio_given = 0; 48 info->fill_ratio_needed = 0; 49 info->factor_mallocs = 0; 50 51 PetscFunctionReturn(0); 52 } 53 54 #undef __FUNC__ 55 #define __FUNC__ "MatScale_SeqDense" 56 int MatScale_SeqDense(Scalar *alpha,Mat inA) 57 { 58 Mat_SeqDense *a = (Mat_SeqDense *) inA->data; 59 int one = 1, nz; 60 61 PetscFunctionBegin; 62 nz = a->m*a->n; 63 BLscal_( &nz, alpha, a->v, &one ); 64 PLogFlops(nz); 65 PetscFunctionReturn(0); 66 } 67 68 /* ---------------------------------------------------------------*/ 69 /* COMMENT: I have chosen to hide column permutation in the pivots, 70 rather than put it in the Mat->col slot.*/ 71 #undef __FUNC__ 72 #define __FUNC__ "MatLUFactor_SeqDense" 73 int MatLUFactor_SeqDense(Mat A,IS row,IS col,double f) 74 { 75 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 76 int info; 77 78 PetscFunctionBegin; 79 if (!mat->pivots) { 80 mat->pivots = (int *) PetscMalloc((mat->m+1)*sizeof(int));CHKPTRQ(mat->pivots); 81 PLogObjectMemory(A,mat->m*sizeof(int)); 82 } 83 LAgetrf_(&mat->m,&mat->n,mat->v,&mat->m,mat->pivots,&info); 84 if (info<0) SETERRQ(PETSC_ERR_LIB,0,"Bad argument to LU factorization"); 85 if (info>0) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,0,"Bad LU factorization"); 86 A->factor = FACTOR_LU; 87 PLogFlops((2*mat->n*mat->n*mat->n)/3); 88 PetscFunctionReturn(0); 89 } 90 91 #undef __FUNC__ 92 #define __FUNC__ "MatConvertSameType_SeqDense" 93 int MatConvertSameType_SeqDense(Mat A,Mat *newmat,int cpvalues) 94 { 95 Mat_SeqDense *mat = (Mat_SeqDense *) A->data, *l; 96 int ierr; 97 Mat newi; 98 99 PetscFunctionBegin; 100 ierr = MatCreateSeqDense(A->comm,mat->m,mat->n,PETSC_NULL,&newi); CHKERRQ(ierr); 101 l = (Mat_SeqDense *) newi->data; 102 if (cpvalues == COPY_VALUES) { 103 PetscMemcpy(l->v,mat->v,mat->m*mat->n*sizeof(Scalar)); 104 } 105 newi->assembled = PETSC_TRUE; 106 *newmat = newi; 107 PetscFunctionReturn(0); 108 } 109 110 #undef __FUNC__ 111 #define __FUNC__ "MatLUFactorSymbolic_SeqDense" 112 int MatLUFactorSymbolic_SeqDense(Mat A,IS row,IS col,double f,Mat *fact) 113 { 114 int ierr; 115 116 PetscFunctionBegin; 117 ierr = MatConvertSameType_SeqDense(A,fact,PETSC_FALSE);CHKERRQ(ierr); 118 PetscFunctionReturn(0); 119 } 120 121 #undef __FUNC__ 122 #define __FUNC__ "MatLUFactorNumeric_SeqDense" 123 int MatLUFactorNumeric_SeqDense(Mat A,Mat *fact) 124 { 125 Mat_SeqDense *mat = (Mat_SeqDense*) A->data, *l = (Mat_SeqDense*) (*fact)->data; 126 int ierr; 127 128 PetscFunctionBegin; 129 /* copy the numerical values */ 130 PetscMemcpy(l->v,mat->v,mat->m*mat->n*sizeof(Scalar)); 131 (*fact)->factor = 0; 132 ierr = MatLUFactor(*fact,0,0,1.0);CHKERRQ(ierr); 133 PetscFunctionReturn(0); 134 } 135 136 #undef __FUNC__ 137 #define __FUNC__ "MatCholeskyFactorSymbolic_SeqDense" 138 int MatCholeskyFactorSymbolic_SeqDense(Mat A,IS row,double f,Mat *fact) 139 { 140 int ierr; 141 142 PetscFunctionBegin; 143 ierr = MatConvert(A,MATSAME,fact);CHKERRQ(ierr); 144 PetscFunctionReturn(0); 145 } 146 147 #undef __FUNC__ 148 #define __FUNC__ "MatCholeskyFactorNumeric_SeqDense" 149 int MatCholeskyFactorNumeric_SeqDense(Mat A,Mat *fact) 150 { 151 int ierr; 152 153 PetscFunctionBegin; 154 ierr = MatCholeskyFactor(*fact,0,1.0);CHKERRQ(ierr); 155 PetscFunctionReturn(0); 156 } 157 158 #undef __FUNC__ 159 #define __FUNC__ "MatCholeskyFactor_SeqDense" 160 int MatCholeskyFactor_SeqDense(Mat A,IS perm,double f) 161 { 162 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 163 int info; 164 165 PetscFunctionBegin; 166 if (mat->pivots) { 167 PetscFree(mat->pivots); 168 PLogObjectMemory(A,-mat->m*sizeof(int)); 169 mat->pivots = 0; 170 } 171 LApotrf_("L",&mat->n,mat->v,&mat->m,&info); 172 if (info) SETERRQ(PETSC_ERR_MAT_CH_ZRPVT,0,"Bad factorization"); 173 A->factor = FACTOR_CHOLESKY; 174 PLogFlops((mat->n*mat->n*mat->n)/3); 175 PetscFunctionReturn(0); 176 } 177 178 #undef __FUNC__ 179 #define __FUNC__ "MatSolve_SeqDense" 180 int MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 181 { 182 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 183 int one = 1, info, ierr; 184 Scalar *x, *y; 185 186 PetscFunctionBegin; 187 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 188 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 189 if (A->factor == FACTOR_LU) { 190 LAgetrs_( "N", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 191 } 192 else if (A->factor == FACTOR_CHOLESKY){ 193 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 194 } 195 else SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Matrix must be factored to solve"); 196 if (info) SETERRQ(PETSC_ERR_LIB,0,"MBad solve"); 197 PLogFlops(mat->n*mat->n - mat->n); 198 PetscFunctionReturn(0); 199 } 200 201 #undef __FUNC__ 202 #define __FUNC__ "MatSolveTrans_SeqDense" 203 int MatSolveTrans_SeqDense(Mat A,Vec xx,Vec yy) 204 { 205 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 206 int one = 1, info; 207 Scalar *x, *y; 208 209 PetscFunctionBegin; 210 VecGetArray(xx,&x); VecGetArray(yy,&y); 211 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 212 /* assume if pivots exist then use LU; else Cholesky */ 213 if (mat->pivots) { 214 LAgetrs_( "T", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 215 } 216 else { 217 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 218 } 219 if (info) SETERRQ(PETSC_ERR_LIB,0,"Bad solve"); 220 PLogFlops(mat->n*mat->n - mat->n); 221 PetscFunctionReturn(0); 222 } 223 224 #undef __FUNC__ 225 #define __FUNC__ "MatSolveAdd_SeqDense" 226 int MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 227 { 228 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 229 int one = 1, info,ierr; 230 Scalar *x, *y, sone = 1.0; 231 Vec tmp = 0; 232 233 PetscFunctionBegin; 234 VecGetArray(xx,&x); VecGetArray(yy,&y); 235 if (yy == zz) { 236 ierr = VecDuplicate(yy,&tmp); CHKERRQ(ierr); 237 PLogObjectParent(A,tmp); 238 ierr = VecCopy(yy,tmp); CHKERRQ(ierr); 239 } 240 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 241 /* assume if pivots exist then use LU; else Cholesky */ 242 if (mat->pivots) { 243 LAgetrs_( "N", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 244 } else { 245 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 246 } 247 if (info) SETERRQ(PETSC_ERR_LIB,0,"Bad solve"); 248 if (tmp) {ierr = VecAXPY(&sone,tmp,yy); CHKERRQ(ierr); ierr = VecDestroy(tmp);CHKERRQ(ierr);} 249 else {ierr = VecAXPY(&sone,zz,yy); CHKERRQ(ierr);} 250 PLogFlops(mat->n*mat->n - mat->n); 251 PetscFunctionReturn(0); 252 } 253 254 #undef __FUNC__ 255 #define __FUNC__ "MatSolveTransAdd_SeqDense" 256 int MatSolveTransAdd_SeqDense(Mat A,Vec xx,Vec zz, Vec yy) 257 { 258 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 259 int one = 1, info,ierr; 260 Scalar *x, *y, sone = 1.0; 261 Vec tmp; 262 263 PetscFunctionBegin; 264 VecGetArray(xx,&x); VecGetArray(yy,&y); 265 if (yy == zz) { 266 ierr = VecDuplicate(yy,&tmp); CHKERRQ(ierr); 267 PLogObjectParent(A,tmp); 268 ierr = VecCopy(yy,tmp); CHKERRQ(ierr); 269 } 270 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 271 /* assume if pivots exist then use LU; else Cholesky */ 272 if (mat->pivots) { 273 LAgetrs_( "T", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 274 } else { 275 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 276 } 277 if (info) SETERRQ(PETSC_ERR_LIB,0,"Bad solve"); 278 if (tmp) { 279 ierr = VecAXPY(&sone,tmp,yy); CHKERRQ(ierr); 280 ierr = VecDestroy(tmp); CHKERRQ(ierr); 281 } else { 282 ierr = VecAXPY(&sone,zz,yy); CHKERRQ(ierr); 283 } 284 PLogFlops(mat->n*mat->n - mat->n); 285 PetscFunctionReturn(0); 286 } 287 /* ------------------------------------------------------------------*/ 288 #undef __FUNC__ 289 #define __FUNC__ "MatRelax_SeqDense" 290 int MatRelax_SeqDense(Mat A,Vec bb,double omega,MatSORType flag, 291 double shift,int its,Vec xx) 292 { 293 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 294 Scalar *x, *b, *v = mat->v, zero = 0.0, xt; 295 int ierr, m = mat->m, i; 296 #if !defined(USE_PETSC_COMPLEX) 297 int o = 1; 298 #endif 299 300 PetscFunctionBegin; 301 if (flag & SOR_ZERO_INITIAL_GUESS) { 302 /* this is a hack fix, should have another version without the second BLdot */ 303 ierr = VecSet(&zero,xx); CHKERRQ(ierr); 304 } 305 VecGetArray(xx,&x); VecGetArray(bb,&b); 306 while (its--) { 307 if (flag & SOR_FORWARD_SWEEP){ 308 for ( i=0; i<m; i++ ) { 309 #if defined(USE_PETSC_COMPLEX) 310 /* cannot use BLAS dot for complex because compiler/linker is 311 not happy about returning a double complex */ 312 int _i; 313 Scalar sum = b[i]; 314 for ( _i=0; _i<m; _i++ ) { 315 sum -= conj(v[i+_i*m])*x[_i]; 316 } 317 xt = sum; 318 #else 319 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 320 #endif 321 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 322 } 323 } 324 if (flag & SOR_BACKWARD_SWEEP) { 325 for ( i=m-1; i>=0; i-- ) { 326 #if defined(USE_PETSC_COMPLEX) 327 /* cannot use BLAS dot for complex because compiler/linker is 328 not happy about returning a double complex */ 329 int _i; 330 Scalar sum = b[i]; 331 for ( _i=0; _i<m; _i++ ) { 332 sum -= conj(v[i+_i*m])*x[_i]; 333 } 334 xt = sum; 335 #else 336 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 337 #endif 338 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 339 } 340 } 341 } 342 PetscFunctionReturn(0); 343 } 344 345 /* -----------------------------------------------------------------*/ 346 #undef __FUNC__ 347 #define __FUNC__ "MatMultTrans_SeqDense" 348 int MatMultTrans_SeqDense(Mat A,Vec xx,Vec yy) 349 { 350 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 351 Scalar *v = mat->v, *x, *y; 352 int _One=1;Scalar _DOne=1.0, _DZero=0.0; 353 354 PetscFunctionBegin; 355 VecGetArray(xx,&x), VecGetArray(yy,&y); 356 LAgemv_("T",&(mat->m),&(mat->n),&_DOne,v,&(mat->m),x,&_One,&_DZero,y,&_One); 357 PLogFlops(2*mat->m*mat->n - mat->n); 358 PetscFunctionReturn(0); 359 } 360 361 #undef __FUNC__ 362 #define __FUNC__ "MatMult_SeqDense" 363 int MatMult_SeqDense(Mat A,Vec xx,Vec yy) 364 { 365 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 366 Scalar *v = mat->v, *x, *y; 367 int _One=1;Scalar _DOne=1.0, _DZero=0.0; 368 369 PetscFunctionBegin; 370 VecGetArray(xx,&x); VecGetArray(yy,&y); 371 LAgemv_( "N", &(mat->m), &(mat->n), &_DOne, v, &(mat->m),x,&_One,&_DZero,y,&_One); 372 PLogFlops(2*mat->m*mat->n - mat->m); 373 PetscFunctionReturn(0); 374 } 375 376 #undef __FUNC__ 377 #define __FUNC__ "MatMultAdd_SeqDense" 378 int MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 379 { 380 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 381 Scalar *v = mat->v, *x, *y, *z; 382 int _One=1; Scalar _DOne=1.0; 383 384 PetscFunctionBegin; 385 VecGetArray(xx,&x); VecGetArray(yy,&y); VecGetArray(zz,&z); 386 if (zz != yy) PetscMemcpy(y,z,mat->m*sizeof(Scalar)); 387 LAgemv_( "N", &(mat->m), &(mat->n),&_DOne,v,&(mat->m),x,&_One,&_DOne,y,&_One); 388 PLogFlops(2*mat->m*mat->n); 389 PetscFunctionReturn(0); 390 } 391 392 #undef __FUNC__ 393 #define __FUNC__ "MatMultTransAdd_SeqDense" 394 int MatMultTransAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 395 { 396 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 397 Scalar *v = mat->v, *x, *y, *z; 398 int _One=1; 399 Scalar _DOne=1.0; 400 401 PetscFunctionBegin; 402 VecGetArray(xx,&x); VecGetArray(yy,&y); VecGetArray(zz,&z); 403 if (zz != yy) PetscMemcpy(y,z,mat->n*sizeof(Scalar)); 404 LAgemv_( "T", &(mat->m), &(mat->n), &_DOne, v, &(mat->m),x,&_One,&_DOne,y,&_One); 405 PLogFlops(2*mat->m*mat->n); 406 PetscFunctionReturn(0); 407 } 408 409 /* -----------------------------------------------------------------*/ 410 #undef __FUNC__ 411 #define __FUNC__ "MatGetRow_SeqDense" 412 int MatGetRow_SeqDense(Mat A,int row,int *ncols,int **cols,Scalar **vals) 413 { 414 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 415 Scalar *v; 416 int i; 417 418 PetscFunctionBegin; 419 *ncols = mat->n; 420 if (cols) { 421 *cols = (int *) PetscMalloc((mat->n+1)*sizeof(int)); CHKPTRQ(*cols); 422 for ( i=0; i<mat->n; i++ ) (*cols)[i] = i; 423 } 424 if (vals) { 425 *vals = (Scalar *) PetscMalloc((mat->n+1)*sizeof(Scalar)); CHKPTRQ(*vals); 426 v = mat->v + row; 427 for ( i=0; i<mat->n; i++ ) {(*vals)[i] = *v; v += mat->m;} 428 } 429 PetscFunctionReturn(0); 430 } 431 432 #undef __FUNC__ 433 #define __FUNC__ "MatRestoreRow_SeqDense" 434 int MatRestoreRow_SeqDense(Mat A,int row,int *ncols,int **cols,Scalar **vals) 435 { 436 if (cols) { PetscFree(*cols); } 437 if (vals) { PetscFree(*vals); } 438 PetscFunctionReturn(0); 439 } 440 /* ----------------------------------------------------------------*/ 441 #undef __FUNC__ 442 #define __FUNC__ "MatSetValues_SeqDense" 443 int MatSetValues_SeqDense(Mat A,int m,int *indexm,int n, 444 int *indexn,Scalar *v,InsertMode addv) 445 { 446 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 447 int i,j; 448 449 PetscFunctionBegin; 450 if (!mat->roworiented) { 451 if (addv == INSERT_VALUES) { 452 for ( j=0; j<n; j++ ) { 453 #if defined(USE_PETSC_BOPT_g) 454 if (indexn[j] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 455 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 456 #endif 457 for ( i=0; i<m; i++ ) { 458 #if defined(USE_PETSC_BOPT_g) 459 if (indexm[i] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 460 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 461 #endif 462 mat->v[indexn[j]*mat->m + indexm[i]] = *v++; 463 } 464 } 465 } else { 466 for ( j=0; j<n; j++ ) { 467 #if defined(USE_PETSC_BOPT_g) 468 if (indexn[j] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 469 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 470 #endif 471 for ( i=0; i<m; i++ ) { 472 #if defined(USE_PETSC_BOPT_g) 473 if (indexm[i] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 474 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 475 #endif 476 mat->v[indexn[j]*mat->m + indexm[i]] += *v++; 477 } 478 } 479 } 480 } else { 481 if (addv == INSERT_VALUES) { 482 for ( i=0; i<m; i++ ) { 483 #if defined(USE_PETSC_BOPT_g) 484 if (indexm[i] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 485 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 486 #endif 487 for ( j=0; j<n; j++ ) { 488 #if defined(USE_PETSC_BOPT_g) 489 if (indexn[j] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 490 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 491 #endif 492 mat->v[indexn[j]*mat->m + indexm[i]] = *v++; 493 } 494 } 495 } else { 496 for ( i=0; i<m; i++ ) { 497 #if defined(USE_PETSC_BOPT_g) 498 if (indexm[i] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 499 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 500 #endif 501 for ( j=0; j<n; j++ ) { 502 #if defined(USE_PETSC_BOPT_g) 503 if (indexn[j] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 504 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 505 #endif 506 mat->v[indexn[j]*mat->m + indexm[i]] += *v++; 507 } 508 } 509 } 510 } 511 PetscFunctionReturn(0); 512 } 513 514 #undef __FUNC__ 515 #define __FUNC__ "MatGetValues_SeqDense" 516 int MatGetValues_SeqDense(Mat A,int m,int *indexm,int n,int *indexn,Scalar *v) 517 { 518 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 519 int i, j; 520 Scalar *vpt = v; 521 522 PetscFunctionBegin; 523 /* row-oriented output */ 524 for ( i=0; i<m; i++ ) { 525 for ( j=0; j<n; j++ ) { 526 *vpt++ = mat->v[indexn[j]*mat->m + indexm[i]]; 527 } 528 } 529 PetscFunctionReturn(0); 530 } 531 532 /* -----------------------------------------------------------------*/ 533 534 #include "sys.h" 535 536 #undef __FUNC__ 537 #define __FUNC__ "MatLoad_SeqDense" 538 int MatLoad_SeqDense(Viewer viewer,MatType type,Mat *A) 539 { 540 Mat_SeqDense *a; 541 Mat B; 542 int *scols, i, j, nz, ierr, fd, header[4], size; 543 int *rowlengths = 0, M, N, *cols; 544 Scalar *vals, *svals, *v; 545 MPI_Comm comm = ((PetscObject)viewer)->comm; 546 547 PetscFunctionBegin; 548 MPI_Comm_size(comm,&size); 549 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"view must have one processor"); 550 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 551 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 552 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"Not matrix object"); 553 M = header[1]; N = header[2]; nz = header[3]; 554 555 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 556 ierr = MatCreateSeqDense(comm,M,N,PETSC_NULL,A); CHKERRQ(ierr); 557 B = *A; 558 a = (Mat_SeqDense *) B->data; 559 560 /* read in nonzero values */ 561 ierr = PetscBinaryRead(fd,a->v,M*N,PETSC_SCALAR); CHKERRQ(ierr); 562 563 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 564 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 565 /* ierr = MatTranspose(B,PETSC_NULL); CHKERRQ(ierr); */ 566 } else { 567 /* read row lengths */ 568 rowlengths = (int*) PetscMalloc( (M+1)*sizeof(int) ); CHKPTRQ(rowlengths); 569 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 570 571 /* create our matrix */ 572 ierr = MatCreateSeqDense(comm,M,N,PETSC_NULL,A); CHKERRQ(ierr); 573 B = *A; 574 a = (Mat_SeqDense *) B->data; 575 v = a->v; 576 577 /* read column indices and nonzeros */ 578 cols = scols = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(cols); 579 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT); CHKERRQ(ierr); 580 vals = svals = (Scalar *) PetscMalloc( (nz+1)*sizeof(Scalar) ); CHKPTRQ(vals); 581 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR); CHKERRQ(ierr); 582 583 /* insert into matrix */ 584 for ( i=0; i<M; i++ ) { 585 for ( j=0; j<rowlengths[i]; j++ ) v[i+M*scols[j]] = svals[j]; 586 svals += rowlengths[i]; scols += rowlengths[i]; 587 } 588 PetscFree(vals); PetscFree(cols); PetscFree(rowlengths); 589 590 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 591 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 592 } 593 PetscFunctionReturn(0); 594 } 595 596 #include "pinclude/pviewer.h" 597 #include "sys.h" 598 599 #undef __FUNC__ 600 #define __FUNC__ "MatView_SeqDense_ASCII" 601 static int MatView_SeqDense_ASCII(Mat A,Viewer viewer) 602 { 603 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 604 int ierr, i, j, format; 605 FILE *fd; 606 char *outputname; 607 Scalar *v; 608 609 PetscFunctionBegin; 610 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 611 ierr = ViewerFileGetOutputname_Private(viewer,&outputname); CHKERRQ(ierr); 612 ierr = ViewerGetFormat(viewer,&format); 613 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 614 PetscFunctionReturn(0); /* do nothing for now */ 615 } 616 else if (format == VIEWER_FORMAT_ASCII_COMMON) { 617 for ( i=0; i<a->m; i++ ) { 618 v = a->v + i; 619 fprintf(fd,"row %d:",i); 620 for ( j=0; j<a->n; j++ ) { 621 #if defined(USE_PETSC_COMPLEX) 622 if (real(*v) != 0.0 && imag(*v) != 0.0) fprintf(fd," %d %g + %g i",j,real(*v),imag(*v)); 623 else if (real(*v)) fprintf(fd," %d %g ",j,real(*v)); 624 #else 625 if (*v) fprintf(fd," %d %g ",j,*v); 626 #endif 627 v += a->m; 628 } 629 fprintf(fd,"\n"); 630 } 631 } else { 632 #if defined(USE_PETSC_COMPLEX) 633 int allreal = 1; 634 /* determine if matrix has all real values */ 635 v = a->v; 636 for ( i=0; i<a->m*a->n; i++ ) { 637 if (imag(v[i])) { allreal = 0; break ;} 638 } 639 #endif 640 for ( i=0; i<a->m; i++ ) { 641 v = a->v + i; 642 for ( j=0; j<a->n; j++ ) { 643 #if defined(USE_PETSC_COMPLEX) 644 if (allreal) { 645 fprintf(fd,"%6.4e ",real(*v)); v += a->m; 646 } else { 647 fprintf(fd,"%6.4e + %6.4e i ",real(*v),imag(*v)); v += a->m; 648 } 649 #else 650 fprintf(fd,"%6.4e ",*v); v += a->m; 651 #endif 652 } 653 fprintf(fd,"\n"); 654 } 655 } 656 fflush(fd); 657 PetscFunctionReturn(0); 658 } 659 660 #undef __FUNC__ 661 #define __FUNC__ "MatView_SeqDense_Binary" 662 static int MatView_SeqDense_Binary(Mat A,Viewer viewer) 663 { 664 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 665 int ict, j, n = a->n, m = a->m, i, fd, *col_lens, ierr, nz = m*n; 666 int format; 667 Scalar *v, *anonz,*vals; 668 669 PetscFunctionBegin; 670 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 671 672 ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr); 673 if (format == VIEWER_FORMAT_BINARY_NATIVE) { 674 /* store the matrix as a dense matrix */ 675 col_lens = (int *) PetscMalloc( 4*sizeof(int) ); CHKPTRQ(col_lens); 676 col_lens[0] = MAT_COOKIE; 677 col_lens[1] = m; 678 col_lens[2] = n; 679 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 680 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,1); CHKERRQ(ierr); 681 PetscFree(col_lens); 682 683 /* write out matrix, by rows */ 684 vals = (Scalar *) PetscMalloc((m*n+1)*sizeof(Scalar)); CHKPTRQ(vals); 685 v = a->v; 686 for ( i=0; i<m; i++ ) { 687 for ( j=0; j<n; j++ ) { 688 vals[i + j*m] = *v++; 689 } 690 } 691 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,0); CHKERRQ(ierr); 692 PetscFree(vals); 693 } else { 694 col_lens = (int *) PetscMalloc( (4+nz)*sizeof(int) ); CHKPTRQ(col_lens); 695 col_lens[0] = MAT_COOKIE; 696 col_lens[1] = m; 697 col_lens[2] = n; 698 col_lens[3] = nz; 699 700 /* store lengths of each row and write (including header) to file */ 701 for ( i=0; i<m; i++ ) col_lens[4+i] = n; 702 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,1); CHKERRQ(ierr); 703 704 /* Possibly should write in smaller increments, not whole matrix at once? */ 705 /* store column indices (zero start index) */ 706 ict = 0; 707 for ( i=0; i<m; i++ ) { 708 for ( j=0; j<n; j++ ) col_lens[ict++] = j; 709 } 710 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,0); CHKERRQ(ierr); 711 PetscFree(col_lens); 712 713 /* store nonzero values */ 714 anonz = (Scalar *) PetscMalloc((nz+1)*sizeof(Scalar)); CHKPTRQ(anonz); 715 ict = 0; 716 for ( i=0; i<m; i++ ) { 717 v = a->v + i; 718 for ( j=0; j<n; j++ ) { 719 anonz[ict++] = *v; v += a->m; 720 } 721 } 722 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,0); CHKERRQ(ierr); 723 PetscFree(anonz); 724 } 725 PetscFunctionReturn(0); 726 } 727 728 #undef __FUNC__ 729 #define __FUNC__ "MatView_SeqDense" 730 int MatView_SeqDense(PetscObject obj,Viewer viewer) 731 { 732 Mat A = (Mat) obj; 733 Mat_SeqDense *a = (Mat_SeqDense*) A->data; 734 ViewerType vtype; 735 int ierr; 736 737 PetscFunctionBegin; 738 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 739 740 if (vtype == MATLAB_VIEWER) { 741 ierr = ViewerMatlabPutArray_Private(viewer,a->m,a->n,a->v); CHKERRQ(ierr); 742 } else if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 743 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 744 } else if (vtype == BINARY_FILE_VIEWER) { 745 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 746 } 747 PetscFunctionReturn(0); 748 } 749 750 #undef __FUNC__ 751 #define __FUNC__ "MatDestroy_SeqDense" 752 int MatDestroy_SeqDense(PetscObject obj) 753 { 754 Mat mat = (Mat) obj; 755 Mat_SeqDense *l = (Mat_SeqDense *) mat->data; 756 int ierr; 757 758 PetscFunctionBegin; 759 #if defined(USE_PETSC_LOG) 760 PLogObjectState(obj,"Rows %d Cols %d",l->m,l->n); 761 #endif 762 if (l->pivots) PetscFree(l->pivots); 763 if (!l->user_alloc) PetscFree(l->v); 764 PetscFree(l); 765 if (mat->mapping) { 766 ierr = ISLocalToGlobalMappingDestroy(mat->mapping); CHKERRQ(ierr); 767 } 768 PLogObjectDestroy(mat); 769 PetscHeaderDestroy(mat); 770 PetscFunctionReturn(0); 771 } 772 773 #undef __FUNC__ 774 #define __FUNC__ "MatTranspose_SeqDense" 775 int MatTranspose_SeqDense(Mat A,Mat *matout) 776 { 777 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 778 int k, j, m, n; 779 Scalar *v, tmp; 780 781 PetscFunctionBegin; 782 v = mat->v; m = mat->m; n = mat->n; 783 if (matout == PETSC_NULL) { /* in place transpose */ 784 if (m != n) { /* malloc temp to hold transpose */ 785 Scalar *w = (Scalar *) PetscMalloc((m+1)*(n+1)*sizeof(Scalar));CHKPTRQ(w); 786 for ( j=0; j<m; j++ ) { 787 for ( k=0; k<n; k++ ) { 788 w[k + j*n] = v[j + k*m]; 789 } 790 } 791 PetscMemcpy(v,w,m*n*sizeof(Scalar)); 792 PetscFree(w); 793 } else { 794 for ( j=0; j<m; j++ ) { 795 for ( k=0; k<j; k++ ) { 796 tmp = v[j + k*n]; 797 v[j + k*n] = v[k + j*n]; 798 v[k + j*n] = tmp; 799 } 800 } 801 } 802 } else { /* out-of-place transpose */ 803 int ierr; 804 Mat tmat; 805 Mat_SeqDense *tmatd; 806 Scalar *v2; 807 ierr = MatCreateSeqDense(A->comm,mat->n,mat->m,PETSC_NULL,&tmat); CHKERRQ(ierr); 808 tmatd = (Mat_SeqDense *) tmat->data; 809 v = mat->v; v2 = tmatd->v; 810 for ( j=0; j<n; j++ ) { 811 for ( k=0; k<m; k++ ) v2[j + k*n] = v[k + j*m]; 812 } 813 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 814 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 815 *matout = tmat; 816 } 817 PetscFunctionReturn(0); 818 } 819 820 #undef __FUNC__ 821 #define __FUNC__ "MatEqual_SeqDense" 822 int MatEqual_SeqDense(Mat A1,Mat A2, PetscTruth *flg) 823 { 824 Mat_SeqDense *mat1 = (Mat_SeqDense *) A1->data; 825 Mat_SeqDense *mat2 = (Mat_SeqDense *) A2->data; 826 int i; 827 Scalar *v1 = mat1->v, *v2 = mat2->v; 828 829 PetscFunctionBegin; 830 if (A2->type != MATSEQDENSE) SETERRQ(PETSC_ERR_SUP,0,"Matrices should be of type MATSEQDENSE"); 831 if (mat1->m != mat2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 832 if (mat1->n != mat2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 833 for ( i=0; i<mat1->m*mat1->n; i++ ) { 834 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 835 v1++; v2++; 836 } 837 *flg = PETSC_TRUE; 838 PetscFunctionReturn(0); 839 } 840 841 #undef __FUNC__ 842 #define __FUNC__ "MatGetDiagonal_SeqDense" 843 int MatGetDiagonal_SeqDense(Mat A,Vec v) 844 { 845 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 846 int i, n, len; 847 Scalar *x, zero = 0.0; 848 849 PetscFunctionBegin; 850 VecSet(&zero,v); 851 VecGetArray(v,&x); VecGetSize(v,&n); 852 len = PetscMin(mat->m,mat->n); 853 if (n != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming mat and vec"); 854 for ( i=0; i<len; i++ ) { 855 x[i] = mat->v[i*mat->m + i]; 856 } 857 PetscFunctionReturn(0); 858 } 859 860 #undef __FUNC__ 861 #define __FUNC__ "MatDiagonalScale_SeqDense" 862 int MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 863 { 864 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 865 Scalar *l,*r,x,*v; 866 int i,j,m = mat->m, n = mat->n; 867 868 PetscFunctionBegin; 869 if (ll) { 870 VecGetArray(ll,&l); VecGetSize(ll,&m); 871 if (m != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vec wrong size"); 872 for ( i=0; i<m; i++ ) { 873 x = l[i]; 874 v = mat->v + i; 875 for ( j=0; j<n; j++ ) { (*v) *= x; v+= m;} 876 } 877 PLogFlops(n*m); 878 } 879 if (rr) { 880 VecGetArray(rr,&r); VecGetSize(rr,&n); 881 if (n != mat->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vec wrong size"); 882 for ( i=0; i<n; i++ ) { 883 x = r[i]; 884 v = mat->v + i*m; 885 for ( j=0; j<m; j++ ) { (*v++) *= x;} 886 } 887 PLogFlops(n*m); 888 } 889 PetscFunctionReturn(0); 890 } 891 892 #undef __FUNC__ 893 #define __FUNC__ "MatNorm_SeqDense" 894 int MatNorm_SeqDense(Mat A,NormType type,double *norm) 895 { 896 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 897 Scalar *v = mat->v; 898 double sum = 0.0; 899 int i, j; 900 901 PetscFunctionBegin; 902 if (type == NORM_FROBENIUS) { 903 for (i=0; i<mat->n*mat->m; i++ ) { 904 #if defined(USE_PETSC_COMPLEX) 905 sum += real(conj(*v)*(*v)); v++; 906 #else 907 sum += (*v)*(*v); v++; 908 #endif 909 } 910 *norm = sqrt(sum); 911 PLogFlops(2*mat->n*mat->m); 912 } else if (type == NORM_1) { 913 *norm = 0.0; 914 for ( j=0; j<mat->n; j++ ) { 915 sum = 0.0; 916 for ( i=0; i<mat->m; i++ ) { 917 sum += PetscAbsScalar(*v); v++; 918 } 919 if (sum > *norm) *norm = sum; 920 } 921 PLogFlops(mat->n*mat->m); 922 } else if (type == NORM_INFINITY) { 923 *norm = 0.0; 924 for ( j=0; j<mat->m; j++ ) { 925 v = mat->v + j; 926 sum = 0.0; 927 for ( i=0; i<mat->n; i++ ) { 928 sum += PetscAbsScalar(*v); v += mat->m; 929 } 930 if (sum > *norm) *norm = sum; 931 } 932 PLogFlops(mat->n*mat->m); 933 } else { 934 SETERRQ(PETSC_ERR_SUP,0,"No two norm"); 935 } 936 PetscFunctionReturn(0); 937 } 938 939 #undef __FUNC__ 940 #define __FUNC__ "MatSetOption_SeqDense" 941 int MatSetOption_SeqDense(Mat A,MatOption op) 942 { 943 Mat_SeqDense *aij = (Mat_SeqDense *) A->data; 944 945 PetscFunctionBegin; 946 if (op == MAT_ROW_ORIENTED) aij->roworiented = 1; 947 else if (op == MAT_COLUMN_ORIENTED) aij->roworiented = 0; 948 else if (op == MAT_ROWS_SORTED || 949 op == MAT_ROWS_UNSORTED || 950 op == MAT_COLUMNS_SORTED || 951 op == MAT_COLUMNS_UNSORTED || 952 op == MAT_SYMMETRIC || 953 op == MAT_STRUCTURALLY_SYMMETRIC || 954 op == MAT_NO_NEW_NONZERO_LOCATIONS || 955 op == MAT_YES_NEW_NONZERO_LOCATIONS || 956 op == MAT_NEW_NONZERO_LOCATION_ERROR || 957 op == MAT_NO_NEW_DIAGONALS || 958 op == MAT_YES_NEW_DIAGONALS || 959 op == MAT_IGNORE_OFF_PROC_ENTRIES || 960 op == MAT_USE_HASH_TABLE) 961 PLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 962 else { 963 SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 964 } 965 PetscFunctionReturn(0); 966 } 967 968 #undef __FUNC__ 969 #define __FUNC__ "MatZeroEntries_SeqDense" 970 int MatZeroEntries_SeqDense(Mat A) 971 { 972 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 973 974 PetscFunctionBegin; 975 PetscMemzero(l->v,l->m*l->n*sizeof(Scalar)); 976 PetscFunctionReturn(0); 977 } 978 979 #undef __FUNC__ 980 #define __FUNC__ "MatGetBlockSize_SeqDense" 981 int MatGetBlockSize_SeqDense(Mat A,int *bs) 982 { 983 PetscFunctionBegin; 984 *bs = 1; 985 PetscFunctionReturn(0); 986 } 987 988 #undef __FUNC__ 989 #define __FUNC__ "MatZeroRows_SeqDense" 990 int MatZeroRows_SeqDense(Mat A,IS is,Scalar *diag) 991 { 992 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 993 int n = l->n, i, j,ierr,N, *rows; 994 Scalar *slot; 995 996 PetscFunctionBegin; 997 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 998 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 999 for ( i=0; i<N; i++ ) { 1000 slot = l->v + rows[i]; 1001 for ( j=0; j<n; j++ ) { *slot = 0.0; slot += n;} 1002 } 1003 if (diag) { 1004 for ( i=0; i<N; i++ ) { 1005 slot = l->v + (n+1)*rows[i]; 1006 *slot = *diag; 1007 } 1008 } 1009 ISRestoreIndices(is,&rows); 1010 PetscFunctionReturn(0); 1011 } 1012 1013 #undef __FUNC__ 1014 #define __FUNC__ "MatGetSize_SeqDense" 1015 int MatGetSize_SeqDense(Mat A,int *m,int *n) 1016 { 1017 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1018 1019 PetscFunctionBegin; 1020 *m = mat->m; *n = mat->n; 1021 PetscFunctionReturn(0); 1022 } 1023 1024 #undef __FUNC__ 1025 #define __FUNC__ "MatGetOwnershipRange_SeqDense" 1026 int MatGetOwnershipRange_SeqDense(Mat A,int *m,int *n) 1027 { 1028 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1029 1030 PetscFunctionBegin; 1031 *m = 0; *n = mat->m; 1032 PetscFunctionReturn(0); 1033 } 1034 1035 #undef __FUNC__ 1036 #define __FUNC__ "MatGetArray_SeqDense" 1037 int MatGetArray_SeqDense(Mat A,Scalar **array) 1038 { 1039 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1040 1041 PetscFunctionBegin; 1042 *array = mat->v; 1043 PetscFunctionReturn(0); 1044 } 1045 1046 #undef __FUNC__ 1047 #define __FUNC__ "MatRestoreArray_SeqDense" 1048 int MatRestoreArray_SeqDense(Mat A,Scalar **array) 1049 { 1050 PetscFunctionBegin; 1051 PetscFunctionReturn(0); 1052 } 1053 1054 #undef __FUNC__ 1055 #define __FUNC__ "MatGetSubMatrix_SeqDense" 1056 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatGetSubMatrixCall scall,Mat *submat) 1057 { 1058 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1059 int nznew, *smap, i, j, ierr, oldcols = mat->n; 1060 int *irow, *icol, nrows, ncols, *cwork; 1061 Scalar *vwork, *val; 1062 Mat newmat; 1063 1064 PetscFunctionBegin; 1065 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1066 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1067 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1068 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1069 1070 smap = (int *) PetscMalloc(oldcols*sizeof(int)); CHKPTRQ(smap); 1071 cwork = (int *) PetscMalloc(ncols*sizeof(int)); CHKPTRQ(cwork); 1072 vwork = (Scalar *) PetscMalloc(ncols*sizeof(Scalar)); CHKPTRQ(vwork); 1073 PetscMemzero((char*)smap,oldcols*sizeof(int)); 1074 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1075 1076 /* Create and fill new matrix */ 1077 ierr = MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat); CHKERRQ(ierr); 1078 for (i=0; i<nrows; i++) { 1079 nznew = 0; 1080 val = mat->v + irow[i]; 1081 for (j=0; j<oldcols; j++) { 1082 if (smap[j]) { 1083 cwork[nznew] = smap[j] - 1; 1084 vwork[nznew++] = val[j * mat->m]; 1085 } 1086 } 1087 ierr = MatSetValues(newmat,1,&i,nznew,cwork,vwork,INSERT_VALUES);CHKERRQ(ierr); 1088 } 1089 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1090 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1091 1092 /* Free work space */ 1093 PetscFree(smap); PetscFree(cwork); PetscFree(vwork); 1094 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1095 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1096 *submat = newmat; 1097 PetscFunctionReturn(0); 1098 } 1099 1100 #undef __FUNC__ 1101 #define __FUNC__ "MatGetSubMatrices_SeqDense" 1102 int MatGetSubMatrices_SeqDense(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall,Mat **B) 1103 { 1104 int ierr,i; 1105 1106 PetscFunctionBegin; 1107 if (scall == MAT_INITIAL_MATRIX) { 1108 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1109 } 1110 1111 for ( i=0; i<n; i++ ) { 1112 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1113 } 1114 PetscFunctionReturn(0); 1115 } 1116 1117 #undef __FUNC__ 1118 #define __FUNC__ "MatCopy_SeqDense" 1119 int MatCopy_SeqDense(Mat A, Mat B) 1120 { 1121 Mat_SeqDense *a = (Mat_SeqDense *) A->data, *b = (Mat_SeqDense *)B->data; 1122 int ierr; 1123 1124 PetscFunctionBegin; 1125 if (B->type != MATSEQDENSE) { 1126 ierr = MatCopy_Basic(A,B);CHKERRQ(ierr); 1127 PetscFunctionReturn(0); 1128 } 1129 if (a->m != b->m || a->n != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"size(B) != size(A)"); 1130 PetscMemcpy(b->v,a->v,a->m*a->n*sizeof(Scalar)); 1131 PetscFunctionReturn(0); 1132 } 1133 1134 /* -------------------------------------------------------------------*/ 1135 static struct _MatOps MatOps = {MatSetValues_SeqDense, 1136 MatGetRow_SeqDense, 1137 MatRestoreRow_SeqDense, 1138 MatMult_SeqDense, 1139 MatMultAdd_SeqDense, 1140 MatMultTrans_SeqDense, 1141 MatMultTransAdd_SeqDense, 1142 MatSolve_SeqDense, 1143 MatSolveAdd_SeqDense, 1144 MatSolveTrans_SeqDense, 1145 MatSolveTransAdd_SeqDense, 1146 MatLUFactor_SeqDense, 1147 MatCholeskyFactor_SeqDense, 1148 MatRelax_SeqDense, 1149 MatTranspose_SeqDense, 1150 MatGetInfo_SeqDense, 1151 MatEqual_SeqDense, 1152 MatGetDiagonal_SeqDense, 1153 MatDiagonalScale_SeqDense, 1154 MatNorm_SeqDense, 1155 0, 1156 0, 1157 0, 1158 MatSetOption_SeqDense, 1159 MatZeroEntries_SeqDense, 1160 MatZeroRows_SeqDense, 1161 MatLUFactorSymbolic_SeqDense, 1162 MatLUFactorNumeric_SeqDense, 1163 MatCholeskyFactorSymbolic_SeqDense, 1164 MatCholeskyFactorNumeric_SeqDense, 1165 MatGetSize_SeqDense, 1166 MatGetSize_SeqDense, 1167 MatGetOwnershipRange_SeqDense, 1168 0, 1169 0, 1170 MatGetArray_SeqDense, 1171 MatRestoreArray_SeqDense, 1172 MatConvertSameType_SeqDense,0,0,0,0, 1173 MatAXPY_SeqDense,MatGetSubMatrices_SeqDense,0, 1174 MatGetValues_SeqDense, 1175 MatCopy_SeqDense,0,MatScale_SeqDense, 1176 0,0,0,MatGetBlockSize_SeqDense}; 1177 1178 #undef __FUNC__ 1179 #define __FUNC__ "MatCreateSeqDense" 1180 /*@C 1181 MatCreateSeqDense - Creates a sequential dense matrix that 1182 is stored in column major order (the usual Fortran 77 manner). Many 1183 of the matrix operations use the BLAS and LAPACK routines. 1184 1185 Input Parameters: 1186 . comm - MPI communicator, set to PETSC_COMM_SELF 1187 . m - number of rows 1188 . n - number of columns 1189 . data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1190 to control all matrix memory allocation. 1191 1192 Output Parameter: 1193 . A - the matrix 1194 1195 Notes: 1196 The data input variable is intended primarily for Fortran programmers 1197 who wish to allocate their own matrix memory space. Most users should 1198 set data=PETSC_NULL. 1199 1200 .keywords: dense, matrix, LAPACK, BLAS 1201 1202 .seealso: MatCreate(), MatSetValues() 1203 @*/ 1204 int MatCreateSeqDense(MPI_Comm comm,int m,int n,Scalar *data,Mat *A) 1205 { 1206 Mat B; 1207 Mat_SeqDense *b; 1208 int ierr,flg,size; 1209 1210 PetscFunctionBegin; 1211 MPI_Comm_size(comm,&size); 1212 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1213 1214 *A = 0; 1215 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQDENSE,comm,MatDestroy,MatView); 1216 PLogObjectCreate(B); 1217 b = (Mat_SeqDense *) PetscMalloc(sizeof(Mat_SeqDense)); CHKPTRQ(b); 1218 PetscMemzero(b,sizeof(Mat_SeqDense)); 1219 PetscMemcpy(B->ops,&MatOps,sizeof(struct _MatOps)); 1220 B->destroy = MatDestroy_SeqDense; 1221 B->view = MatView_SeqDense; 1222 B->factor = 0; 1223 B->mapping = 0; 1224 PLogObjectMemory(B,sizeof(struct _p_Mat)); 1225 B->data = (void *) b; 1226 1227 b->m = m; B->m = m; B->M = m; 1228 b->n = n; B->n = n; B->N = n; 1229 b->pivots = 0; 1230 b->roworiented = 1; 1231 if (data == PETSC_NULL) { 1232 b->v = (Scalar*) PetscMalloc((m*n+1)*sizeof(Scalar)); CHKPTRQ(b->v); 1233 PetscMemzero(b->v,m*n*sizeof(Scalar)); 1234 b->user_alloc = 0; 1235 PLogObjectMemory(B,n*m*sizeof(Scalar)); 1236 } 1237 else { /* user-allocated storage */ 1238 b->v = data; 1239 b->user_alloc = 1; 1240 } 1241 ierr = OptionsHasName(PETSC_NULL,"-help",&flg); CHKERRQ(ierr); 1242 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr);} 1243 1244 *A = B; 1245 PetscFunctionReturn(0); 1246 } 1247 1248 1249 1250 1251 1252 1253