1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: dense.c,v 1.141 1998/04/03 22:01:04 balay 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 } 784 PetscFunctionReturn(0); 785 } 786 787 #undef __FUNC__ 788 #define __FUNC__ "MatDestroy_SeqDense" 789 int MatDestroy_SeqDense(Mat mat) 790 { 791 Mat_SeqDense *l = (Mat_SeqDense *) mat->data; 792 int ierr; 793 794 PetscFunctionBegin; 795 #if defined(USE_PETSC_LOG) 796 PLogObjectState((PetscObject)mat,"Rows %d Cols %d",l->m,l->n); 797 #endif 798 if (l->pivots) PetscFree(l->pivots); 799 if (!l->user_alloc) PetscFree(l->v); 800 PetscFree(l); 801 if (mat->mapping) { 802 ierr = ISLocalToGlobalMappingDestroy(mat->mapping); CHKERRQ(ierr); 803 } 804 PLogObjectDestroy(mat); 805 PetscHeaderDestroy(mat); 806 PetscFunctionReturn(0); 807 } 808 809 #undef __FUNC__ 810 #define __FUNC__ "MatTranspose_SeqDense" 811 int MatTranspose_SeqDense(Mat A,Mat *matout) 812 { 813 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 814 int k, j, m, n; 815 Scalar *v, tmp; 816 817 PetscFunctionBegin; 818 v = mat->v; m = mat->m; n = mat->n; 819 if (matout == PETSC_NULL) { /* in place transpose */ 820 if (m != n) { /* malloc temp to hold transpose */ 821 Scalar *w = (Scalar *) PetscMalloc((m+1)*(n+1)*sizeof(Scalar));CHKPTRQ(w); 822 for ( j=0; j<m; j++ ) { 823 for ( k=0; k<n; k++ ) { 824 w[k + j*n] = v[j + k*m]; 825 } 826 } 827 PetscMemcpy(v,w,m*n*sizeof(Scalar)); 828 PetscFree(w); 829 } else { 830 for ( j=0; j<m; j++ ) { 831 for ( k=0; k<j; k++ ) { 832 tmp = v[j + k*n]; 833 v[j + k*n] = v[k + j*n]; 834 v[k + j*n] = tmp; 835 } 836 } 837 } 838 } else { /* out-of-place transpose */ 839 int ierr; 840 Mat tmat; 841 Mat_SeqDense *tmatd; 842 Scalar *v2; 843 ierr = MatCreateSeqDense(A->comm,mat->n,mat->m,PETSC_NULL,&tmat); CHKERRQ(ierr); 844 tmatd = (Mat_SeqDense *) tmat->data; 845 v = mat->v; v2 = tmatd->v; 846 for ( j=0; j<n; j++ ) { 847 for ( k=0; k<m; k++ ) v2[j + k*n] = v[k + j*m]; 848 } 849 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 850 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 851 *matout = tmat; 852 } 853 PetscFunctionReturn(0); 854 } 855 856 #undef __FUNC__ 857 #define __FUNC__ "MatEqual_SeqDense" 858 int MatEqual_SeqDense(Mat A1,Mat A2, PetscTruth *flg) 859 { 860 Mat_SeqDense *mat1 = (Mat_SeqDense *) A1->data; 861 Mat_SeqDense *mat2 = (Mat_SeqDense *) A2->data; 862 int i; 863 Scalar *v1 = mat1->v, *v2 = mat2->v; 864 865 PetscFunctionBegin; 866 if (A2->type != MATSEQDENSE) SETERRQ(PETSC_ERR_SUP,0,"Matrices should be of type MATSEQDENSE"); 867 if (mat1->m != mat2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 868 if (mat1->n != mat2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 869 for ( i=0; i<mat1->m*mat1->n; i++ ) { 870 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 871 v1++; v2++; 872 } 873 *flg = PETSC_TRUE; 874 PetscFunctionReturn(0); 875 } 876 877 #undef __FUNC__ 878 #define __FUNC__ "MatGetDiagonal_SeqDense" 879 int MatGetDiagonal_SeqDense(Mat A,Vec v) 880 { 881 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 882 int ierr,i, n, len; 883 Scalar *x, zero = 0.0; 884 885 PetscFunctionBegin; 886 ierr = VecSet(&zero,v);CHKERRQ(ierr); 887 ierr = VecGetArray(v,&x); CHKERRQ(ierr); 888 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 889 len = PetscMin(mat->m,mat->n); 890 if (n != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming mat and vec"); 891 for ( i=0; i<len; i++ ) { 892 x[i] = mat->v[i*mat->m + i]; 893 } 894 ierr = VecRestoreArray(v,&x); CHKERRQ(ierr); 895 PetscFunctionReturn(0); 896 } 897 898 #undef __FUNC__ 899 #define __FUNC__ "MatDiagonalScale_SeqDense" 900 int MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 901 { 902 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 903 Scalar *l,*r,x,*v; 904 int ierr,i,j,m = mat->m, n = mat->n; 905 906 PetscFunctionBegin; 907 if (ll) { 908 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 909 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 910 if (m != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vec wrong size"); 911 for ( i=0; i<m; i++ ) { 912 x = l[i]; 913 v = mat->v + i; 914 for ( j=0; j<n; j++ ) { (*v) *= x; v+= m;} 915 } 916 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 917 PLogFlops(n*m); 918 } 919 if (rr) { 920 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 921 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 922 if (n != mat->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vec wrong size"); 923 for ( i=0; i<n; i++ ) { 924 x = r[i]; 925 v = mat->v + i*m; 926 for ( j=0; j<m; j++ ) { (*v++) *= x;} 927 } 928 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 929 PLogFlops(n*m); 930 } 931 PetscFunctionReturn(0); 932 } 933 934 #undef __FUNC__ 935 #define __FUNC__ "MatNorm_SeqDense" 936 int MatNorm_SeqDense(Mat A,NormType type,double *norm) 937 { 938 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 939 Scalar *v = mat->v; 940 double sum = 0.0; 941 int i, j; 942 943 PetscFunctionBegin; 944 if (type == NORM_FROBENIUS) { 945 for (i=0; i<mat->n*mat->m; i++ ) { 946 #if defined(USE_PETSC_COMPLEX) 947 sum += real(conj(*v)*(*v)); v++; 948 #else 949 sum += (*v)*(*v); v++; 950 #endif 951 } 952 *norm = sqrt(sum); 953 PLogFlops(2*mat->n*mat->m); 954 } else if (type == NORM_1) { 955 *norm = 0.0; 956 for ( j=0; j<mat->n; j++ ) { 957 sum = 0.0; 958 for ( i=0; i<mat->m; i++ ) { 959 sum += PetscAbsScalar(*v); v++; 960 } 961 if (sum > *norm) *norm = sum; 962 } 963 PLogFlops(mat->n*mat->m); 964 } else if (type == NORM_INFINITY) { 965 *norm = 0.0; 966 for ( j=0; j<mat->m; j++ ) { 967 v = mat->v + j; 968 sum = 0.0; 969 for ( i=0; i<mat->n; i++ ) { 970 sum += PetscAbsScalar(*v); v += mat->m; 971 } 972 if (sum > *norm) *norm = sum; 973 } 974 PLogFlops(mat->n*mat->m); 975 } else { 976 SETERRQ(PETSC_ERR_SUP,0,"No two norm"); 977 } 978 PetscFunctionReturn(0); 979 } 980 981 #undef __FUNC__ 982 #define __FUNC__ "MatSetOption_SeqDense" 983 int MatSetOption_SeqDense(Mat A,MatOption op) 984 { 985 Mat_SeqDense *aij = (Mat_SeqDense *) A->data; 986 987 PetscFunctionBegin; 988 if (op == MAT_ROW_ORIENTED) aij->roworiented = 1; 989 else if (op == MAT_COLUMN_ORIENTED) aij->roworiented = 0; 990 else if (op == MAT_ROWS_SORTED || 991 op == MAT_ROWS_UNSORTED || 992 op == MAT_COLUMNS_SORTED || 993 op == MAT_COLUMNS_UNSORTED || 994 op == MAT_SYMMETRIC || 995 op == MAT_STRUCTURALLY_SYMMETRIC || 996 op == MAT_NO_NEW_NONZERO_LOCATIONS || 997 op == MAT_YES_NEW_NONZERO_LOCATIONS || 998 op == MAT_NEW_NONZERO_LOCATION_ERROR || 999 op == MAT_NO_NEW_DIAGONALS || 1000 op == MAT_YES_NEW_DIAGONALS || 1001 op == MAT_IGNORE_OFF_PROC_ENTRIES || 1002 op == MAT_USE_HASH_TABLE) 1003 PLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1004 else { 1005 SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 1006 } 1007 PetscFunctionReturn(0); 1008 } 1009 1010 #undef __FUNC__ 1011 #define __FUNC__ "MatZeroEntries_SeqDense" 1012 int MatZeroEntries_SeqDense(Mat A) 1013 { 1014 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 1015 1016 PetscFunctionBegin; 1017 PetscMemzero(l->v,l->m*l->n*sizeof(Scalar)); 1018 PetscFunctionReturn(0); 1019 } 1020 1021 #undef __FUNC__ 1022 #define __FUNC__ "MatGetBlockSize_SeqDense" 1023 int MatGetBlockSize_SeqDense(Mat A,int *bs) 1024 { 1025 PetscFunctionBegin; 1026 *bs = 1; 1027 PetscFunctionReturn(0); 1028 } 1029 1030 #undef __FUNC__ 1031 #define __FUNC__ "MatZeroRows_SeqDense" 1032 int MatZeroRows_SeqDense(Mat A,IS is,Scalar *diag) 1033 { 1034 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 1035 int n = l->n, i, j,ierr,N, *rows; 1036 Scalar *slot; 1037 1038 PetscFunctionBegin; 1039 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 1040 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 1041 for ( i=0; i<N; i++ ) { 1042 slot = l->v + rows[i]; 1043 for ( j=0; j<n; j++ ) { *slot = 0.0; slot += n;} 1044 } 1045 if (diag) { 1046 for ( i=0; i<N; i++ ) { 1047 slot = l->v + (n+1)*rows[i]; 1048 *slot = *diag; 1049 } 1050 } 1051 ISRestoreIndices(is,&rows); 1052 PetscFunctionReturn(0); 1053 } 1054 1055 #undef __FUNC__ 1056 #define __FUNC__ "MatGetSize_SeqDense" 1057 int MatGetSize_SeqDense(Mat A,int *m,int *n) 1058 { 1059 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1060 1061 PetscFunctionBegin; 1062 *m = mat->m; *n = mat->n; 1063 PetscFunctionReturn(0); 1064 } 1065 1066 #undef __FUNC__ 1067 #define __FUNC__ "MatGetOwnershipRange_SeqDense" 1068 int MatGetOwnershipRange_SeqDense(Mat A,int *m,int *n) 1069 { 1070 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1071 1072 PetscFunctionBegin; 1073 *m = 0; *n = mat->m; 1074 PetscFunctionReturn(0); 1075 } 1076 1077 #undef __FUNC__ 1078 #define __FUNC__ "MatGetArray_SeqDense" 1079 int MatGetArray_SeqDense(Mat A,Scalar **array) 1080 { 1081 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1082 1083 PetscFunctionBegin; 1084 *array = mat->v; 1085 PetscFunctionReturn(0); 1086 } 1087 1088 #undef __FUNC__ 1089 #define __FUNC__ "MatRestoreArray_SeqDense" 1090 int MatRestoreArray_SeqDense(Mat A,Scalar **array) 1091 { 1092 PetscFunctionBegin; 1093 PetscFunctionReturn(0); 1094 } 1095 1096 #undef __FUNC__ 1097 #define __FUNC__ "MatGetSubMatrix_SeqDense" 1098 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatGetSubMatrixCall scall,Mat *submat) 1099 { 1100 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1101 int nznew, *smap, i, j, ierr, oldcols = mat->n; 1102 int *irow, *icol, nrows, ncols, *cwork; 1103 Scalar *vwork, *val; 1104 Mat newmat; 1105 1106 PetscFunctionBegin; 1107 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1108 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1109 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1110 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1111 1112 smap = (int *) PetscMalloc(oldcols*sizeof(int)); CHKPTRQ(smap); 1113 cwork = (int *) PetscMalloc(ncols*sizeof(int)); CHKPTRQ(cwork); 1114 vwork = (Scalar *) PetscMalloc(ncols*sizeof(Scalar)); CHKPTRQ(vwork); 1115 PetscMemzero((char*)smap,oldcols*sizeof(int)); 1116 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1117 1118 /* Create and fill new matrix */ 1119 ierr = MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat); CHKERRQ(ierr); 1120 for (i=0; i<nrows; i++) { 1121 nznew = 0; 1122 val = mat->v + irow[i]; 1123 for (j=0; j<oldcols; j++) { 1124 if (smap[j]) { 1125 cwork[nznew] = smap[j] - 1; 1126 vwork[nznew++] = val[j * mat->m]; 1127 } 1128 } 1129 ierr = MatSetValues(newmat,1,&i,nznew,cwork,vwork,INSERT_VALUES);CHKERRQ(ierr); 1130 } 1131 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1132 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1133 1134 /* Free work space */ 1135 PetscFree(smap); PetscFree(cwork); PetscFree(vwork); 1136 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1137 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1138 *submat = newmat; 1139 PetscFunctionReturn(0); 1140 } 1141 1142 #undef __FUNC__ 1143 #define __FUNC__ "MatGetSubMatrices_SeqDense" 1144 int MatGetSubMatrices_SeqDense(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall,Mat **B) 1145 { 1146 int ierr,i; 1147 1148 PetscFunctionBegin; 1149 if (scall == MAT_INITIAL_MATRIX) { 1150 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1151 } 1152 1153 for ( i=0; i<n; i++ ) { 1154 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1155 } 1156 PetscFunctionReturn(0); 1157 } 1158 1159 #undef __FUNC__ 1160 #define __FUNC__ "MatCopy_SeqDense" 1161 int MatCopy_SeqDense(Mat A, Mat B) 1162 { 1163 Mat_SeqDense *a = (Mat_SeqDense *) A->data, *b = (Mat_SeqDense *)B->data; 1164 int ierr; 1165 1166 PetscFunctionBegin; 1167 if (B->type != MATSEQDENSE) { 1168 ierr = MatCopy_Basic(A,B);CHKERRQ(ierr); 1169 PetscFunctionReturn(0); 1170 } 1171 if (a->m != b->m || a->n != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"size(B) != size(A)"); 1172 PetscMemcpy(b->v,a->v,a->m*a->n*sizeof(Scalar)); 1173 PetscFunctionReturn(0); 1174 } 1175 1176 /* -------------------------------------------------------------------*/ 1177 static struct _MatOps MatOps = {MatSetValues_SeqDense, 1178 MatGetRow_SeqDense, 1179 MatRestoreRow_SeqDense, 1180 MatMult_SeqDense, 1181 MatMultAdd_SeqDense, 1182 MatMultTrans_SeqDense, 1183 MatMultTransAdd_SeqDense, 1184 MatSolve_SeqDense, 1185 MatSolveAdd_SeqDense, 1186 MatSolveTrans_SeqDense, 1187 MatSolveTransAdd_SeqDense, 1188 MatLUFactor_SeqDense, 1189 MatCholeskyFactor_SeqDense, 1190 MatRelax_SeqDense, 1191 MatTranspose_SeqDense, 1192 MatGetInfo_SeqDense, 1193 MatEqual_SeqDense, 1194 MatGetDiagonal_SeqDense, 1195 MatDiagonalScale_SeqDense, 1196 MatNorm_SeqDense, 1197 0, 1198 0, 1199 0, 1200 MatSetOption_SeqDense, 1201 MatZeroEntries_SeqDense, 1202 MatZeroRows_SeqDense, 1203 MatLUFactorSymbolic_SeqDense, 1204 MatLUFactorNumeric_SeqDense, 1205 MatCholeskyFactorSymbolic_SeqDense, 1206 MatCholeskyFactorNumeric_SeqDense, 1207 MatGetSize_SeqDense, 1208 MatGetSize_SeqDense, 1209 MatGetOwnershipRange_SeqDense, 1210 0, 1211 0, 1212 MatGetArray_SeqDense, 1213 MatRestoreArray_SeqDense, 1214 MatConvertSameType_SeqDense,0,0,0,0, 1215 MatAXPY_SeqDense,MatGetSubMatrices_SeqDense,0, 1216 MatGetValues_SeqDense, 1217 MatCopy_SeqDense,0,MatScale_SeqDense, 1218 0,0,0,MatGetBlockSize_SeqDense}; 1219 1220 #undef __FUNC__ 1221 #define __FUNC__ "MatCreateSeqDense" 1222 /*@C 1223 MatCreateSeqDense - Creates a sequential dense matrix that 1224 is stored in column major order (the usual Fortran 77 manner). Many 1225 of the matrix operations use the BLAS and LAPACK routines. 1226 1227 Input Parameters: 1228 . comm - MPI communicator, set to PETSC_COMM_SELF 1229 . m - number of rows 1230 . n - number of columns 1231 . data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1232 to control all matrix memory allocation. 1233 1234 Output Parameter: 1235 . A - the matrix 1236 1237 Notes: 1238 The data input variable is intended primarily for Fortran programmers 1239 who wish to allocate their own matrix memory space. Most users should 1240 set data=PETSC_NULL. 1241 1242 .keywords: dense, matrix, LAPACK, BLAS 1243 1244 .seealso: MatCreate(), MatSetValues() 1245 @*/ 1246 int MatCreateSeqDense(MPI_Comm comm,int m,int n,Scalar *data,Mat *A) 1247 { 1248 Mat B; 1249 Mat_SeqDense *b; 1250 int ierr,flg,size; 1251 1252 PetscFunctionBegin; 1253 MPI_Comm_size(comm,&size); 1254 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1255 1256 *A = 0; 1257 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQDENSE,comm,MatDestroy,MatView); 1258 PLogObjectCreate(B); 1259 b = (Mat_SeqDense *) PetscMalloc(sizeof(Mat_SeqDense)); CHKPTRQ(b); 1260 PetscMemzero(b,sizeof(Mat_SeqDense)); 1261 PetscMemcpy(B->ops,&MatOps,sizeof(struct _MatOps)); 1262 B->ops->destroy = MatDestroy_SeqDense; 1263 B->ops->view = MatView_SeqDense; 1264 B->factor = 0; 1265 B->mapping = 0; 1266 PLogObjectMemory(B,sizeof(struct _p_Mat)); 1267 B->data = (void *) b; 1268 1269 b->m = m; B->m = m; B->M = m; 1270 b->n = n; B->n = n; B->N = n; 1271 b->pivots = 0; 1272 b->roworiented = 1; 1273 if (data == PETSC_NULL) { 1274 b->v = (Scalar*) PetscMalloc((m*n+1)*sizeof(Scalar)); CHKPTRQ(b->v); 1275 PetscMemzero(b->v,m*n*sizeof(Scalar)); 1276 b->user_alloc = 0; 1277 PLogObjectMemory(B,n*m*sizeof(Scalar)); 1278 } 1279 else { /* user-allocated storage */ 1280 b->v = data; 1281 b->user_alloc = 1; 1282 } 1283 ierr = OptionsHasName(PETSC_NULL,"-help",&flg); CHKERRQ(ierr); 1284 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr);} 1285 1286 *A = B; 1287 PetscFunctionReturn(0); 1288 } 1289 1290 1291 1292 1293 1294 1295