1 2 /* 3 Standard Fortran style matrices 4 */ 5 #include "ptscimpl.h" 6 #include "plapack.h" 7 #include "matimpl.h" 8 #include "math.h" 9 #include "vec/vecimpl.h" 10 11 typedef struct { 12 Scalar *v; 13 int roworiented; 14 int m,n,pad; 15 int *pivots; /* pivots in LU factorization */ 16 } MatiSD; 17 18 19 static int MatiSDnz(Mat matin,int *nz) 20 { 21 MatiSD *mat = (MatiSD *) matin->data; 22 int i,N = mat->m*mat->n,count = 0; 23 Scalar *v = mat->v; 24 for ( i=0; i<N; i++ ) {if (*v != 0.0) count++; v++;} 25 *nz = count; return 0; 26 } 27 static int MatiSDmemory(Mat matin,int *mem) 28 { 29 MatiSD *mat = (MatiSD *) matin->data; 30 *mem = mat->m*mat->n*sizeof(Scalar); return 0; 31 } 32 33 /* ---------------------------------------------------------------*/ 34 /* COMMENT: I have chosen to hide column permutation in the pivots, 35 rather than put it in the Mat->col slot.*/ 36 static int MatiSDlufactor(Mat matin,IS row,IS col) 37 { 38 MatiSD *mat = (MatiSD *) matin->data; 39 int ierr, one = 1, info; 40 if (!mat->pivots) { 41 mat->pivots = (int *) MALLOC( mat->m*sizeof(int) ); 42 CHKPTR(mat->pivots); 43 } 44 LAgetrf_(&mat->m,&mat->n,mat->v,&mat->m,mat->pivots,&info); 45 if (info) SETERR(1,"Bad LU factorization"); 46 matin->factor = FACTOR_LU; 47 return 0; 48 } 49 static int MatiSDlufactorsymbolic(Mat matin,IS row,IS col,Mat *fact) 50 { 51 int ierr; 52 if (ierr = MatCopy(matin,fact)) SETERR(ierr,0); 53 return 0; 54 } 55 static int MatiSDlufactornumeric(Mat matin,Mat *fact) 56 { 57 return MatLUFactor(*fact,0,0); 58 } 59 static int MatiSDchfactorsymbolic(Mat matin,IS row,Mat *fact) 60 { 61 int ierr; 62 if (ierr = MatCopy(matin,fact)) SETERR(ierr,0); 63 return 0; 64 } 65 static int MatiSDchfactornumeric(Mat matin,Mat *fact) 66 { 67 return MatCholeskyFactor(*fact,0); 68 } 69 static int MatiSDchfactor(Mat matin,IS perm) 70 { 71 MatiSD *mat = (MatiSD *) matin->data; 72 int ierr, one = 1, info; 73 if (mat->pivots) {FREE(mat->pivots); mat->pivots = 0;} 74 LApotrf_("L",&mat->n,mat->v,&mat->m,&info); 75 if (info) SETERR(1,"Bad Cholesky factorization"); 76 matin->factor = FACTOR_CHOLESKY; 77 return 0; 78 } 79 80 static int MatiSDsolve(Mat matin,Vec xx,Vec yy) 81 { 82 MatiSD *mat = (MatiSD *) matin->data; 83 int i,j, one = 1, info; 84 Scalar *v = mat->v, *x, *y; 85 VecGetArray(xx,&x); VecGetArray(yy,&y); 86 MEMCPY(y,x,mat->m*sizeof(Scalar)); 87 /* assume if pivots exist then LU else Cholesky */ 88 if (mat->pivots) { 89 LAgetrs_( "N", &mat->m, &one, mat->v, &mat->m, mat->pivots, 90 y, &mat->m, &info ); 91 } 92 else { 93 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m, 94 y, &mat->m, &info ); 95 } 96 if (info) SETERR(1,"Bad solve"); 97 return 0; 98 } 99 static int MatiSDsolvetrans(Mat matin,Vec xx,Vec yy) 100 { 101 MatiSD *mat = (MatiSD *) matin->data; 102 int i,j, one = 1, info; 103 Scalar *v = mat->v, *x, *y; 104 VecGetArray(xx,&x); VecGetArray(yy,&y); 105 MEMCPY(y,x,mat->m*sizeof(Scalar)); 106 /* assume if pivots exist then LU else Cholesky */ 107 if (mat->pivots) { 108 LAgetrs_( "T", &mat->m, &one, mat->v, &mat->m, mat->pivots, 109 y, &mat->m, &info ); 110 } 111 else { 112 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m, 113 y, &mat->m, &info ); 114 } 115 if (info) SETERR(1,"Bad solve"); 116 return 0; 117 } 118 static int MatiSDsolveadd(Mat matin,Vec xx,Vec zz,Vec yy) 119 { 120 MatiSD *mat = (MatiSD *) matin->data; 121 int i,j, one = 1, info,ierr; 122 Scalar *v = mat->v, *x, *y, sone = 1.0; 123 Vec tmp = 0; 124 VecGetArray(xx,&x); VecGetArray(yy,&y); 125 if (yy == zz) { 126 ierr = VecCreate(yy,&tmp); CHKERR(ierr); 127 ierr = VecCopy(yy,tmp); CHKERR(ierr); 128 } 129 MEMCPY(y,x,mat->m*sizeof(Scalar)); 130 /* assume if pivots exist then LU else Cholesky */ 131 if (mat->pivots) { 132 LAgetrs_( "N", &mat->m, &one, mat->v, &mat->m, mat->pivots, 133 y, &mat->m, &info ); 134 } 135 else { 136 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m, 137 y, &mat->m, &info ); 138 } 139 if (info) SETERR(1,"Bad solve"); 140 if (tmp) {VecAXPY(&sone,tmp,yy); VecDestroy(tmp);} 141 else VecAXPY(&sone,zz,yy); 142 return 0; 143 } 144 static int MatiSDsolvetransadd(Mat matin,Vec xx,Vec zz, Vec yy) 145 { 146 MatiSD *mat = (MatiSD *) matin->data; 147 int i,j, one = 1, info,ierr; 148 Scalar *v = mat->v, *x, *y, sone = 1.0; 149 Vec tmp; 150 VecGetArray(xx,&x); VecGetArray(yy,&y); 151 if (yy == zz) { 152 ierr = VecCreate(yy,&tmp); CHKERR(ierr); 153 ierr = VecCopy(yy,tmp); CHKERR(ierr); 154 } 155 MEMCPY(y,x,mat->m*sizeof(Scalar)); 156 /* assume if pivots exist then LU else Cholesky */ 157 if (mat->pivots) { 158 LAgetrs_( "T", &mat->m, &one, mat->v, &mat->m, mat->pivots, 159 y, &mat->m, &info ); 160 } 161 else { 162 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m, 163 y, &mat->m, &info ); 164 } 165 if (info) SETERR(1,"Bad solve"); 166 if (tmp) {VecAXPY(&sone,tmp,yy); VecDestroy(tmp);} 167 else VecAXPY(&sone,zz,yy); 168 return 0; 169 } 170 /* ------------------------------------------------------------------*/ 171 static int MatiSDrelax(Mat matin,Vec bb,double omega,int flag,double shift, 172 int its,Vec xx) 173 { 174 MatiSD *mat = (MatiSD *) matin->data; 175 Scalar *x, *b, *v = mat->v, zero = 0.0, xt; 176 int o = 1, tmp,n = mat->n,ierr, m = mat->m, i, j; 177 178 if (flag & SOR_ZERO_INITIAL_GUESS) { 179 /* this is a hack fix, should have another version without 180 the second BLdot */ 181 if (ierr = VecSet(&zero,xx)) SETERR(ierr,0); 182 } 183 VecGetArray(xx,&x); VecGetArray(bb,&b); 184 while (its--) { 185 if (flag & SOR_FORWARD_SWEEP){ 186 for ( i=0; i<m; i++ ) { 187 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 188 x[i] = (1. - omega)*x[i] + omega*(xt/(v[i + i*m]+shift) + x[i]); 189 } 190 } 191 if (flag & SOR_BACKWARD_SWEEP) { 192 for ( i=m-1; i>=0; i-- ) { 193 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 194 x[i] = (1. - omega)*x[i] + omega*(xt/(v[i + i*m]+shift) + x[i]); 195 } 196 } 197 } 198 return 0; 199 } 200 201 /* -----------------------------------------------------------------*/ 202 static int MatiSDmulttrans(Mat matin,Vec xx,Vec yy) 203 { 204 MatiSD *mat = (MatiSD *) matin->data; 205 Scalar *v = mat->v, *x, *y; 206 int _One=1;Scalar _DOne=1.0, _DZero=0.0; 207 VecGetArray(xx,&x), VecGetArray(yy,&y); 208 LAgemv_( "T", &(mat->m), &(mat->n), &_DOne, v, &(mat->m), 209 x, &_One, &_DZero, y, &_One ); 210 return 0; 211 } 212 static int MatiSDmult(Mat matin,Vec xx,Vec yy) 213 { 214 MatiSD *mat = (MatiSD *) matin->data; 215 Scalar *v = mat->v, *x, *y; 216 int _One=1;Scalar _DOne=1.0, _DZero=0.0; 217 VecGetArray(xx,&x); VecGetArray(yy,&y); 218 LAgemv_( "N", &(mat->m), &(mat->n), &_DOne, v, &(mat->m), 219 x, &_One, &_DZero, y, &_One ); 220 return 0; 221 } 222 static int MatiSDmultadd(Mat matin,Vec xx,Vec zz,Vec yy) 223 { 224 MatiSD *mat = (MatiSD *) matin->data; 225 Scalar *v = mat->v, *x, *y, *z; 226 int _One=1; Scalar _DOne=1.0, _DZero=0.0; 227 VecGetArray(xx,&x); VecGetArray(yy,&y); VecGetArray(zz,&z); 228 if (zz != yy) MEMCPY(y,z,mat->m*sizeof(Scalar)); 229 LAgemv_( "N", &(mat->m), &(mat->n), &_DOne, v, &(mat->m), 230 x, &_One, &_DOne, y, &_One ); 231 return 0; 232 } 233 static int MatiSDmulttransadd(Mat matin,Vec xx,Vec zz,Vec yy) 234 { 235 MatiSD *mat = (MatiSD *) matin->data; 236 Scalar *v = mat->v, *x, *y, *z; 237 int _One=1; 238 Scalar _DOne=1.0, _DZero=0.0; 239 VecGetArray(xx,&x); VecGetArray(yy,&y); 240 VecGetArray(zz,&z); 241 if (zz != yy) MEMCPY(y,z,mat->m*sizeof(Scalar)); 242 LAgemv_( "T", &(mat->m), &(mat->n), &_DOne, v, &(mat->m), 243 x, &_One, &_DOne, y, &_One ); 244 return 0; 245 } 246 247 /* -----------------------------------------------------------------*/ 248 static int MatiSDgetrow(Mat matin,int row,int *ncols,int **cols, 249 Scalar **vals) 250 { 251 MatiSD *mat = (MatiSD *) matin->data; 252 Scalar *v; 253 int i; 254 *ncols = mat->n; 255 if (cols) { 256 *cols = (int *) MALLOC(mat->n*sizeof(int)); CHKPTR(*cols); 257 for ( i=0; i<mat->n; i++ ) *cols[i] = i; 258 } 259 if (vals) { 260 *vals = (Scalar *) MALLOC(mat->n*sizeof(Scalar)); CHKPTR(*vals); 261 v = mat->v + row; 262 for ( i=0; i<mat->n; i++ ) {*vals[i] = *v; v += mat->m;} 263 } 264 return 0; 265 } 266 static int MatiSDrestorerow(Mat matin,int row,int *ncols,int **cols, 267 Scalar **vals) 268 { 269 MatiSD *mat = (MatiSD *) matin->data; 270 if (cols) { FREE(*cols); } 271 if (vals) { FREE(*vals); } 272 return 0; 273 } 274 /* ----------------------------------------------------------------*/ 275 static int MatiSDinsert(Mat matin,int m,int *indexm,int n, 276 int *indexn,Scalar *v,InsertMode addv) 277 { 278 MatiSD *mat = (MatiSD *) matin->data; 279 int i,j; 280 281 if (!mat->roworiented) { 282 if (addv == InsertValues) { 283 for ( j=0; j<n; j++ ) { 284 if (indexn[j] < 0) {*v += m; continue;} 285 for ( i=0; i<m; i++ ) { 286 if (indexm[i] < 0) {*v++; continue;} 287 mat->v[indexn[j]*mat->m + indexm[i]] = *v++; 288 } 289 } 290 } 291 else { 292 for ( j=0; j<n; j++ ) { 293 if (indexn[j] < 0) {*v += m; continue;} 294 for ( i=0; i<m; i++ ) { 295 if (indexm[i] < 0) {*v++; continue;} 296 mat->v[indexn[j]*mat->m + indexm[i]] += *v++; 297 } 298 } 299 } 300 } 301 else { 302 if (addv == InsertValues) { 303 for ( i=0; i<m; i++ ) { 304 if (indexm[i] < 0) {*v += n; continue;} 305 for ( j=0; j<n; j++ ) { 306 if (indexn[j] < 0) {*v++; continue;} 307 mat->v[indexn[j]*mat->m + indexm[i]] = *v++; 308 } 309 } 310 } 311 else { 312 for ( i=0; i<m; i++ ) { 313 if (indexm[i] < 0) {*v += n; continue;} 314 for ( j=0; j<n; j++ ) { 315 if (indexn[j] < 0) {*v++; continue;} 316 mat->v[indexn[j]*mat->m + indexm[i]] += *v++; 317 } 318 } 319 } 320 } 321 return 0; 322 } 323 324 /* -----------------------------------------------------------------*/ 325 static int MatiSDcopy(Mat matin,Mat *newmat) 326 { 327 MatiSD *mat = (MatiSD *) matin->data; 328 int ierr; 329 Mat newi; 330 MatiSD *l; 331 if (ierr = MatCreateSequentialDense(mat->m,mat->n,&newi)) SETERR(ierr,0); 332 l = (MatiSD *) newi->data; 333 MEMCPY(l->v,mat->v,mat->m*mat->n*sizeof(Scalar)); 334 *newmat = newi; 335 return 0; 336 } 337 #include "viewer.h" 338 339 int MatiSDview(PetscObject obj,Viewer ptr) 340 { 341 Mat matin = (Mat) obj; 342 MatiSD *mat = (MatiSD *) matin->data; 343 Scalar *v; 344 int i,j; 345 PetscObject ojb = (PetscObject) ptr; 346 347 if (obj && obj->cookie == VIEWER_COOKIE && obj->type == MATLAB_VIEWER) { 348 return ViewerMatlabPutArray(ptr,mat->m,mat->n,mat->v); 349 } 350 else { 351 for ( i=0; i<mat->m; i++ ) { 352 v = mat->v + i; 353 for ( j=0; j<mat->n; j++ ) { 354 #if defined(PETSC_COMPLEX) 355 printf("%6.4e + %6.4e i ",real(*v),imag(*v)); v += mat->m; 356 #else 357 printf("%6.4e ",*v); v += mat->m; 358 #endif 359 } 360 printf("\n"); 361 } 362 } 363 return 0; 364 } 365 366 367 static int MatiSDdestroy(PetscObject obj) 368 { 369 Mat mat = (Mat) obj; 370 MatiSD *l = (MatiSD *) mat->data; 371 if (l->pivots) FREE(l->pivots); 372 FREE(l); 373 FREE(mat); 374 return 0; 375 } 376 377 static int MatiSDtrans(Mat matin) 378 { 379 MatiSD *mat = (MatiSD *) matin->data; 380 int k,j; 381 Scalar *v = mat->v, tmp; 382 if (mat->m != mat->n) { 383 SETERR(1,"Cannot transpose rectangular dense matrix"); 384 } 385 for ( j=0; j<mat->m; j++ ) { 386 for ( k=0; k<j; k++ ) { 387 tmp = v[j + k*mat->n]; 388 v[j + k*mat->n] = v[k + j*mat->n]; 389 v[k + j*mat->n] = tmp; 390 } 391 } 392 return 0; 393 } 394 395 static int MatiSDequal(Mat matin1,Mat matin2) 396 { 397 MatiSD *mat1 = (MatiSD *) matin1->data; 398 MatiSD *mat2 = (MatiSD *) matin2->data; 399 int i; 400 Scalar *v1 = mat1->v, *v2 = mat2->v; 401 if (mat1->m != mat2->m) return 0; 402 if (mat1->n != mat2->n) return 0; 403 for ( i=0; i<mat1->m*mat1->n; i++ ) { 404 if (*v1 != *v2) return 0; 405 v1++; v2++; 406 } 407 return 1; 408 } 409 410 static int MatiSDgetdiag(Mat matin,Vec v) 411 { 412 MatiSD *mat = (MatiSD *) matin->data; 413 int i,j, n; 414 Scalar *x, zero = 0.0; 415 CHKTYPE(v,SEQVECTOR); 416 VecGetArray(v,&x); VecGetSize(v,&n); 417 if (n != mat->m) SETERR(1,"Nonconforming matrix and vector"); 418 for ( i=0; i<mat->m; i++ ) { 419 x[i] = mat->v[i*mat->m + i]; 420 } 421 return 0; 422 } 423 424 static int MatiSDscale(Mat matin,Vec ll,Vec rr) 425 { 426 MatiSD *mat = (MatiSD *) matin->data; 427 Scalar *l,*r,x,*v; 428 int i,j,m = mat->m, n = mat->n; 429 if (l) { 430 VecGetArray(ll,&l); VecGetSize(ll,&m); 431 if (m != mat->m) SETERR(1,"Left scaling vector wrong length"); 432 for ( i=0; i<m; i++ ) { 433 x = l[i]; 434 v = mat->v + i; 435 for ( j=0; j<n; j++ ) { (*v) *= x; v+= m;} 436 } 437 } 438 if (l) { 439 VecGetArray(rr,&r); VecGetSize(rr,&n); 440 if (n != mat->n) SETERR(1,"Right scaling vector wrong length"); 441 for ( i=0; i<n; i++ ) { 442 x = r[i]; 443 v = mat->v + i*m; 444 for ( j=0; j<m; j++ ) { (*v++) *= x;} 445 } 446 } 447 return 0; 448 } 449 450 451 static int MatiSDnorm(Mat matin,int type,double *norm) 452 { 453 MatiSD *mat = (MatiSD *) matin->data; 454 Scalar *v = mat->v; 455 double sum = 0.0; 456 int i, j; 457 if (type == NORM_FROBENIUS) { 458 for (i=0; i<mat->n*mat->m; i++ ) { 459 #if defined(PETSC_COMPLEX) 460 sum += real(conj(*v)*(*v)); v++; 461 #else 462 sum += (*v)*(*v); v++; 463 #endif 464 } 465 *norm = sqrt(sum); 466 } 467 else if (type == NORM_1) { 468 *norm = 0.0; 469 for ( j=0; j<mat->n; j++ ) { 470 sum = 0.0; 471 for ( i=0; i<mat->m; i++ ) { 472 #if defined(PETSC_COMPLEX) 473 sum += abs(*v++); 474 #else 475 sum += fabs(*v++); 476 #endif 477 } 478 if (sum > *norm) *norm = sum; 479 } 480 } 481 else if (type == NORM_INFINITY) { 482 *norm = 0.0; 483 for ( j=0; j<mat->m; j++ ) { 484 v = mat->v + j; 485 sum = 0.0; 486 for ( i=0; i<mat->n; i++ ) { 487 #if defined(PETSC_COMPLEX) 488 sum += abs(*v); v += mat->m; 489 #else 490 sum += fabs(*v); v += mat->m; 491 #endif 492 } 493 if (sum > *norm) *norm = sum; 494 } 495 } 496 else { 497 SETERR(1,"No support for the two norm yet"); 498 } 499 return 0; 500 } 501 502 static int MatiDenseinsopt(Mat aijin,int op) 503 { 504 MatiSD *aij = (MatiSD *) aijin->data; 505 if (op == ROW_ORIENTED) aij->roworiented = 1; 506 else if (op == COLUMN_ORIENTED) aij->roworiented = 0; 507 /* doesn't care about sorted rows or columns */ 508 return 0; 509 } 510 511 static int MatiZero(Mat A) 512 { 513 MatiSD *l = (MatiSD *) A->data; 514 MEMSET(l->v,0,l->m*l->n*sizeof(Scalar)); 515 return 0; 516 } 517 518 static int MatiZerorows(Mat A,IS is,Scalar *diag) 519 { 520 MatiSD *l = (MatiSD *) A->data; 521 int m = l->m, n = l->n, i, j,ierr,N, *rows; 522 Scalar *slot; 523 ierr = ISGetLocalSize(is,&N); CHKERR(ierr); 524 ierr = ISGetIndices(is,&rows); CHKERR(ierr); 525 for ( i=0; i<N; i++ ) { 526 slot = l->v + rows[i]; 527 for ( j=0; j<n; j++ ) { *slot = 0.0; slot += n;} 528 } 529 if (diag) { 530 for ( i=0; i<N; i++ ) { 531 slot = l->v + (n+1)*rows[i]; 532 *slot = *diag; 533 } 534 } 535 ISRestoreIndices(is,&rows); 536 return 0; 537 } 538 /* -------------------------------------------------------------------*/ 539 static struct _MatOps MatOps = {MatiSDinsert, 540 MatiSDgetrow, MatiSDrestorerow, 541 MatiSDmult, MatiSDmultadd, MatiSDmulttrans, MatiSDmulttransadd, 542 MatiSDsolve,MatiSDsolveadd,MatiSDsolvetrans,MatiSDsolvetransadd, 543 MatiSDlufactor,MatiSDchfactor, 544 MatiSDrelax, 545 MatiSDtrans, 546 MatiSDnz,MatiSDmemory,MatiSDequal, 547 MatiSDcopy, 548 MatiSDgetdiag,MatiSDscale,MatiSDnorm, 549 0,0, 550 0, MatiDenseinsopt,MatiZero,MatiZerorows,0, 551 MatiSDlufactorsymbolic,MatiSDlufactornumeric, 552 MatiSDchfactorsymbolic,MatiSDchfactornumeric 553 }; 554 /*@ 555 MatCreateSequentialDense - Creates a sequential dense matrix that 556 is stored in the usual Fortran 77 manner. Many of the matrix 557 operations use the BLAS and LAPACK routines. 558 559 Input Parameters: 560 . m, n - the number of rows and columns in the matrix. 561 562 Output Parameter: 563 . newmat - the matrix created. 564 565 Keywords: dense matrix, lapack, blas 566 @*/ 567 int MatCreateSequentialDense(int m,int n,Mat *newmat) 568 { 569 int size = sizeof(MatiSD) + m*n*sizeof(Scalar); 570 Mat mat; 571 MatiSD *l; 572 *newmat = 0; 573 CREATEHEADER(mat,_Mat); 574 l = (MatiSD *) MALLOC(size); CHKPTR(l); 575 mat->cookie = MAT_COOKIE; 576 mat->ops = &MatOps; 577 mat->destroy = MatiSDdestroy; 578 mat->view = MatiSDview; 579 mat->data = (void *) l; 580 mat->type = MATDENSESEQ; 581 mat->factor = 0; 582 mat->col = 0; 583 mat->row = 0; 584 585 l->m = m; 586 l->n = n; 587 l->v = (Scalar *) (l + 1); 588 l->pivots = 0; 589 l->roworiented = 1; 590 591 MEMSET(l->v,0,m*n*sizeof(Scalar)); 592 *newmat = mat; 593 return 0; 594 } 595 596 int MatiSDCreate(Mat matin,Mat *newmat) 597 { 598 MatiSD *m = (MatiSD *) matin->data; 599 return MatCreateSequentialDense(m->m,m->n,newmat); 600 } 601