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 {return 0;} 101 102 /* ------------------------------------------------------------------*/ 103 static int MatiSDrelax(Mat matin,Vec bb,double omega,int flag,double shift, 104 int its,Vec xx) 105 { 106 MatiSD *mat = (MatiSD *) matin->data; 107 Scalar *x, *b, *v = mat->v, zero = 0.0, xt; 108 int o = 1, tmp,n = mat->n,ierr, m = mat->m, i, j; 109 110 if (flag & SOR_ZERO_INITIAL_GUESS) { 111 /* this is a hack fix, should have another version without 112 the second BLdot */ 113 if (ierr = VecSet(&zero,xx)) SETERR(ierr,0); 114 } 115 VecGetArray(xx,&x); VecGetArray(bb,&b); 116 while (its--) { 117 if (flag & SOR_FORWARD_SWEEP){ 118 for ( i=0; i<m; i++ ) { 119 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 120 x[i] = (1. - omega)*x[i] + omega*(xt/(v[i + i*m]+shift) + x[i]); 121 } 122 } 123 if (flag & SOR_BACKWARD_SWEEP) { 124 for ( i=m-1; i>=0; i-- ) { 125 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 126 x[i] = (1. - omega)*x[i] + omega*(xt/(v[i + i*m]+shift) + x[i]); 127 } 128 } 129 } 130 return 0; 131 } 132 133 /* -----------------------------------------------------------------*/ 134 static int MatiSDmulttrans(Mat matin,Vec xx,Vec yy) 135 { 136 MatiSD *mat = (MatiSD *) matin->data; 137 Scalar *v = mat->v, *x, *y; 138 int _One=1;Scalar _DOne=1.0, _DZero=0.0; 139 VecGetArray(xx,&x), VecGetArray(yy,&y); 140 LAgemv_( "T", &(mat->m), &(mat->n), &_DOne, v, &(mat->m), 141 x, &_One, &_DZero, y, &_One ); 142 return 0; 143 } 144 static int MatiSDmult(Mat matin,Vec xx,Vec yy) 145 { 146 MatiSD *mat = (MatiSD *) matin->data; 147 Scalar *v = mat->v, *x, *y; 148 int _One=1;Scalar _DOne=1.0, _DZero=0.0; 149 VecGetArray(xx,&x); VecGetArray(yy,&y); 150 LAgemv_( "N", &(mat->m), &(mat->n), &_DOne, v, &(mat->m), 151 x, &_One, &_DZero, y, &_One ); 152 return 0; 153 } 154 static int MatiSDmultadd(Mat matin,Vec xx,Vec zz,Vec yy) 155 { 156 MatiSD *mat = (MatiSD *) matin->data; 157 Scalar *v = mat->v, *x, *y, *z; 158 int _One=1; Scalar _DOne=1.0, _DZero=0.0; 159 VecGetArray(xx,&x); VecGetArray(yy,&y); VecGetArray(zz,&z); 160 if (zz != yy) MEMCPY(y,z,mat->m*sizeof(Scalar)); 161 LAgemv_( "N", &(mat->m), &(mat->n), &_DOne, v, &(mat->m), 162 x, &_One, &_DOne, y, &_One ); 163 return 0; 164 } 165 static int MatiSDmulttransadd(Mat matin,Vec xx,Vec zz,Vec yy) 166 { 167 MatiSD *mat = (MatiSD *) matin->data; 168 Scalar *v = mat->v, *x, *y, *z; 169 int _One=1; 170 Scalar _DOne=1.0, _DZero=0.0; 171 VecGetArray(xx,&x); VecGetArray(yy,&y); 172 VecGetArray(zz,&z); 173 if (zz != yy) MEMCPY(y,z,mat->m*sizeof(Scalar)); 174 LAgemv_( "T", &(mat->m), &(mat->n), &_DOne, v, &(mat->m), 175 x, &_One, &_DOne, y, &_One ); 176 return 0; 177 } 178 179 /* -----------------------------------------------------------------*/ 180 static int MatiSDgetrow(Mat matin,int row,int *ncols,int **cols, 181 Scalar **vals) 182 { 183 MatiSD *mat = (MatiSD *) matin->data; 184 Scalar *v; 185 int i; 186 *ncols = mat->n; 187 if (cols) { 188 *cols = (int *) MALLOC(mat->n*sizeof(int)); CHKPTR(*cols); 189 for ( i=0; i<mat->n; i++ ) *cols[i] = i; 190 } 191 if (vals) { 192 *vals = (Scalar *) MALLOC(mat->n*sizeof(Scalar)); CHKPTR(*vals); 193 v = mat->v + row; 194 for ( i=0; i<mat->n; i++ ) {*vals[i] = *v; v += mat->m;} 195 } 196 return 0; 197 } 198 static int MatiSDrestorerow(Mat matin,int row,int *ncols,int **cols, 199 Scalar **vals) 200 { 201 MatiSD *mat = (MatiSD *) matin->data; 202 if (cols) { FREE(*cols); } 203 if (vals) { FREE(*vals); } 204 return 0; 205 } 206 /* ----------------------------------------------------------------*/ 207 static int MatiSDinsert(Mat matin,int m,int *indexm,int n, 208 int *indexn,Scalar *v,InsertMode addv) 209 { 210 MatiSD *mat = (MatiSD *) matin->data; 211 int i,j; 212 213 if (!mat->roworiented) { 214 if (addv == InsertValues) { 215 for ( j=0; j<n; j++ ) { 216 if (indexn[j] < 0) {*v += m; continue;} 217 for ( i=0; i<m; i++ ) { 218 if (indexm[i] < 0) {*v++; continue;} 219 mat->v[indexn[j]*mat->m + indexm[i]] = *v++; 220 } 221 } 222 } 223 else { 224 for ( j=0; j<n; j++ ) { 225 if (indexn[j] < 0) {*v += m; continue;} 226 for ( i=0; i<m; i++ ) { 227 if (indexm[i] < 0) {*v++; continue;} 228 mat->v[indexn[j]*mat->m + indexm[i]] += *v++; 229 } 230 } 231 } 232 } 233 else { 234 if (addv == InsertValues) { 235 for ( i=0; i<m; i++ ) { 236 if (indexm[i] < 0) {*v += n; continue;} 237 for ( j=0; j<n; j++ ) { 238 if (indexn[j] < 0) {*v++; continue;} 239 mat->v[indexn[j]*mat->m + indexm[i]] = *v++; 240 } 241 } 242 } 243 else { 244 for ( i=0; i<m; i++ ) { 245 if (indexm[i] < 0) {*v += n; continue;} 246 for ( j=0; j<n; j++ ) { 247 if (indexn[j] < 0) {*v++; continue;} 248 mat->v[indexn[j]*mat->m + indexm[i]] += *v++; 249 } 250 } 251 } 252 } 253 return 0; 254 } 255 256 /* -----------------------------------------------------------------*/ 257 static int MatiSDcopy(Mat matin,Mat *newmat) 258 { 259 MatiSD *mat = (MatiSD *) matin->data; 260 int ierr; 261 Mat newi; 262 MatiSD *l; 263 if (ierr = MatCreateSequentialDense(mat->m,mat->n,&newi)) SETERR(ierr,0); 264 l = (MatiSD *) newi->data; 265 MEMCPY(l->v,mat->v,mat->m*mat->n*sizeof(Scalar)); 266 *newmat = newi; 267 return 0; 268 } 269 270 int MatiSDview(PetscObject obj,Viewer ptr) 271 { 272 Mat matin = (Mat) obj; 273 MatiSD *mat = (MatiSD *) matin->data; 274 Scalar *v; 275 int i,j; 276 for ( i=0; i<mat->m; i++ ) { 277 v = mat->v + i; 278 for ( j=0; j<mat->n; j++ ) { 279 #if defined(PETSC_COMPLEX) 280 printf("%6.4e + %6.4e i ",real(*v),imag(*v)); v += mat->m; 281 #else 282 printf("%6.4e ",*v); v += mat->m; 283 #endif 284 } 285 printf("\n"); 286 } 287 return 0; 288 } 289 290 291 static int MatiSDdestroy(PetscObject obj) 292 { 293 Mat mat = (Mat) obj; 294 MatiSD *l = (MatiSD *) mat->data; 295 if (l->pivots) FREE(l->pivots); 296 FREE(l); 297 FREE(mat); 298 return 0; 299 } 300 301 static int MatiSDtrans(Mat matin) 302 { 303 MatiSD *mat = (MatiSD *) matin->data; 304 int k,j; 305 Scalar *v = mat->v, tmp; 306 if (mat->m != mat->n) { 307 SETERR(1,"Cannot transpose rectangular dense matrix"); 308 } 309 for ( j=0; j<mat->m; j++ ) { 310 for ( k=0; k<j; k++ ) { 311 tmp = v[j + k*mat->n]; 312 v[j + k*mat->n] = v[k + j*mat->n]; 313 v[k + j*mat->n] = tmp; 314 } 315 } 316 return 0; 317 } 318 319 static int MatiSDequal(Mat matin1,Mat matin2) 320 { 321 MatiSD *mat1 = (MatiSD *) matin1->data; 322 MatiSD *mat2 = (MatiSD *) matin2->data; 323 int i; 324 Scalar *v1 = mat1->v, *v2 = mat2->v; 325 if (mat1->m != mat2->m) return 0; 326 if (mat1->n != mat2->n) return 0; 327 for ( i=0; i<mat1->m*mat1->n; i++ ) { 328 if (*v1 != *v2) return 0; 329 v1++; v2++; 330 } 331 return 1; 332 } 333 334 static int MatiSDgetdiag(Mat matin,Vec v) 335 { 336 MatiSD *mat = (MatiSD *) matin->data; 337 int i,j, n; 338 Scalar *x, zero = 0.0; 339 CHKTYPE(v,SEQVECTOR); 340 VecGetArray(v,&x); VecGetSize(v,&n); 341 if (n != mat->m) SETERR(1,"Nonconforming matrix and vector"); 342 for ( i=0; i<mat->m; i++ ) { 343 x[i] = mat->v[i*mat->m + i]; 344 } 345 return 0; 346 } 347 348 static int MatiSDscale(Mat matin,Vec l,Vec r) 349 { 350 return 0; 351 } 352 353 static int MatiSDnorm(Mat matin,int type,double *norm) 354 { 355 MatiSD *mat = (MatiSD *) matin->data; 356 Scalar *v = mat->v; 357 double sum = 0.0; 358 int i, j; 359 if (type == NORM_FROBENIUS) { 360 for (i=0; i<mat->n*mat->m; i++ ) { 361 #if defined(PETSC_COMPLEX) 362 sum += real(conj(*v)*(*v)); v++; 363 #else 364 sum += (*v)*(*v); v++; 365 #endif 366 } 367 *norm = sqrt(sum); 368 } 369 else if (type == NORM_1) { 370 *norm = 0.0; 371 for ( j=0; j<mat->n; j++ ) { 372 sum = 0.0; 373 for ( i=0; i<mat->m; i++ ) { 374 #if defined(PETSC_COMPLEX) 375 sum += abs(*v++); 376 #else 377 sum += fabs(*v++); 378 #endif 379 } 380 if (sum > *norm) *norm = sum; 381 } 382 } 383 else if (type == NORM_INFINITY) { 384 *norm = 0.0; 385 for ( j=0; j<mat->m; j++ ) { 386 v = mat->v + j; 387 sum = 0.0; 388 for ( i=0; i<mat->n; i++ ) { 389 #if defined(PETSC_COMPLEX) 390 sum += abs(*v); v += mat->m; 391 #else 392 sum += fabs(*v); v += mat->m; 393 #endif 394 } 395 if (sum > *norm) *norm = sum; 396 } 397 } 398 else { 399 SETERR(1,"No support for the two norm yet"); 400 } 401 return 0; 402 } 403 404 static int MatiDenseinsopt(Mat aijin,int op) 405 { 406 MatiSD *aij = (MatiSD *) aijin->data; 407 if (op == ROW_ORIENTED) aij->roworiented = 1; 408 else if (op == COLUMN_ORIENTED) aij->roworiented = 0; 409 /* doesn't care about sorted rows or columns */ 410 return 0; 411 } 412 413 static int MatiZero(Mat A) 414 { 415 MatiSD *l = (MatiSD *) A->data; 416 MEMSET(l->v,0,l->m*l->n*sizeof(Scalar)); 417 return 0; 418 } 419 420 static int MatiZerorows(Mat A,IS is,Scalar *diag) 421 { 422 MatiSD *l = (MatiSD *) A->data; 423 int m = l->m, n = l->n, i, j,ierr,N, *rows; 424 Scalar *slot; 425 ierr = ISGetLocalSize(is,&N); CHKERR(ierr); 426 ierr = ISGetIndices(is,&rows); CHKERR(ierr); 427 for ( i=0; i<N; i++ ) { 428 slot = l->v + rows[i]; 429 for ( j=0; j<n; j++ ) { *slot = 0.0; slot += n;} 430 } 431 if (diag) { 432 for ( i=0; i<N; i++ ) { 433 slot = l->v + (n+1)*rows[i]; 434 *slot = *diag; 435 } 436 } 437 ISRestoreIndices(is,&rows); 438 return 0; 439 } 440 /* -------------------------------------------------------------------*/ 441 static struct _MatOps MatOps = {MatiSDinsert, 442 MatiSDgetrow, MatiSDrestorerow, 443 MatiSDmult, MatiSDmultadd, MatiSDmulttrans, MatiSDmulttransadd, 444 MatiSDsolve,0,MatiSDsolvetrans,0, 445 MatiSDlufactor,MatiSDchfactor, 446 MatiSDrelax, 447 MatiSDtrans, 448 MatiSDnz,MatiSDmemory,MatiSDequal, 449 MatiSDcopy, 450 MatiSDgetdiag,MatiSDscale,MatiSDnorm, 451 0,0, 452 0, MatiDenseinsopt,MatiZero,MatiZerorows,0, 453 MatiSDlufactorsymbolic,MatiSDlufactornumeric, 454 MatiSDchfactorsymbolic,MatiSDchfactornumeric 455 }; 456 /*@ 457 MatCreateSequentialDense - Creates a sequential dense matrix that 458 is stored in the usual Fortran 77 manner. Many of the matrix 459 operations use the BLAS and LAPACK routines. 460 461 Input Parameters: 462 . m, n - the number of rows and columns in the matrix. 463 464 Output Parameter: 465 . newmat - the matrix created. 466 467 Keywords: dense matrix, lapack, blas 468 @*/ 469 int MatCreateSequentialDense(int m,int n,Mat *newmat) 470 { 471 int size = sizeof(MatiSD) + m*n*sizeof(Scalar); 472 Mat mat; 473 MatiSD *l; 474 *newmat = 0; 475 CREATEHEADER(mat,_Mat); 476 l = (MatiSD *) MALLOC(size); CHKPTR(l); 477 mat->cookie = MAT_COOKIE; 478 mat->ops = &MatOps; 479 mat->destroy = MatiSDdestroy; 480 mat->view = MatiSDview; 481 mat->data = (void *) l; 482 mat->type = MATDENSESEQ; 483 mat->factor = 0; 484 mat->col = 0; 485 mat->row = 0; 486 mat->outofrange= 0; 487 mat->Mlow = 0; 488 mat->Mhigh = m; 489 mat->Nlow = 0; 490 mat->Nhigh = n; 491 492 l->m = m; 493 l->n = n; 494 l->v = (Scalar *) (l + 1); 495 l->pivots = 0; 496 l->roworiented = 1; 497 498 MEMSET(l->v,0,m*n*sizeof(Scalar)); 499 *newmat = mat; 500 return 0; 501 } 502 503 int MatiSDCreate(Mat matin,Mat *newmat) 504 { 505 MatiSD *m = (MatiSD *) matin->data; 506 return MatCreateSequentialDense(m->m,m->n,newmat); 507 } 508