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