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