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