1 2 #include "aij.h" 3 #include "vec/vecimpl.h" 4 #include "inline/spops.h" 5 6 int SpToSymmetricIJ(Matiaij*,int**,int**); 7 int SpOrderND(int,int*,int*,int*); 8 int SpOrder1WD(int,int*,int*,int*); 9 int SpOrderQMD(int,int*,int*,int*); 10 int SpOrderRCM(int,int*,int*,int*); 11 12 static int MatAIJreorder(Mat mat,int type,IS *rperm, IS *cperm) 13 { 14 Matiaij *aij = (Matiaij *) mat->data; 15 int i, ierr, *ia, *ja, *perma; 16 17 perma = (int *) MALLOC( aij->n*sizeof(int) ); CHKPTR(perma); 18 19 if (ierr = SpToSymmetricIJ( aij, &ia, &ja )) SETERR(ierr,0); 20 21 if (type == ORDER_NATURAL) { 22 for ( i=0; i<aij->n; i++ ) perma[i] = i; 23 } 24 else if (type == ORDER_ND) { 25 ierr = SpOrderND( aij->n, ia, ja, perma ); 26 } 27 else if (type == ORDER_1WD) { 28 ierr = SpOrder1WD( aij->n, ia, ja, perma ); 29 } 30 else if (type == ORDER_RCM) { 31 ierr = SpOrderRCM( aij->n, ia, ja, perma ); 32 } 33 else if (type == ORDER_QMD) { 34 ierr = SpOrderQMD( aij->n, ia, ja, perma ); 35 } 36 else SETERR(1,"Cannot performing ordering requested"); 37 if (ierr) SETERR(ierr,0); 38 FREE(ia); FREE(ja); 39 40 if (ierr = ISCreateSequential(aij->n,perma,rperm)) SETERR(ierr,0); 41 ISSetPermutation(*rperm); 42 FREE(perma); 43 *cperm = *rperm; /* so far all permutations are symmetric.*/ 44 return 0; 45 } 46 47 48 #define CHUNCKSIZE 5 49 50 /* This version has row oriented v */ 51 static int MatiAIJAddValues(Mat matin,Scalar *v,int m,int *idxm,int n, 52 int *idxn) 53 { 54 Matiaij *mat = (Matiaij *) matin->data; 55 int *rp,k,a,b,t,ii,row,nrow,i,col,l,rmax, N, sorted = mat->sorted; 56 int *imax = mat->imax, *ai = mat->i, *ailen = mat->ilen; 57 int *aj = mat->j, nonew = mat->nonew; 58 Scalar *ap,value, *aa = mat->a; 59 60 for ( k=0; k<m; k++ ) { /* loop over added rows */ 61 row = idxm[k]; rp = aj + ai[row] - 1; ap = aa + ai[row] - 1; 62 rmax = imax[row]; nrow = ailen[row]; 63 a = 0; 64 for ( l=0; l<n; l++ ) { /* loop over added columns */ 65 col = idxn[l] + 1; value = *v++; 66 if (nrow) { 67 if (!sorted) a = 0; b = nrow; 68 while (b-a > 5) { 69 t = (b+a)/2; 70 if (rp[t] > col) b = t; 71 else a = t; 72 } 73 for ( i=a; i<b; i++ ) { 74 if (rp[i] > col) break; 75 if (rp[i] == col) {ap[i] += value; goto noinsert;} 76 } 77 if (nonew) goto noinsert; 78 if (nrow >= rmax) { 79 /* there is no extra room in row, therefore enlarge */ 80 int new_nz = ai[mat->m] + CHUNCKSIZE,len,*new_i,*new_j; 81 Scalar *new_a; 82 fprintf(stderr,"Warning, enlarging matrix storage \n"); 83 84 /* malloc new storage space */ 85 len = new_nz*(sizeof(int)+sizeof(Scalar))+(mat->m+1)*sizeof(int); 86 new_a = (Scalar *) MALLOC( len ); CHKPTR(new_a); 87 new_j = (int *) (new_a + new_nz); 88 new_i = new_j + new_nz; 89 90 /* copy over old data into new slots */ 91 for ( ii=0; ii<row+1; ii++ ) {new_i[ii] = ai[ii];} 92 for ( ii=row+1; ii<mat->m+1; ii++ ) {new_i[ii] = ai[ii]+CHUNCKSIZE;} 93 MEMCPY(new_j,aj,(ai[row]+nrow-1)*sizeof(int)); 94 len = (new_nz - CHUNCKSIZE - ai[row] - nrow + 1); 95 MEMCPY(new_j+ai[row]-1+nrow+CHUNCKSIZE,aj+ai[row]-1+nrow, 96 len*sizeof(int)); 97 MEMCPY(new_a,aa,(ai[row]+nrow-1)*sizeof(Scalar)); 98 MEMCPY(new_a+ai[row]-1+nrow+CHUNCKSIZE,aa+ai[row]-1+nrow, 99 len*sizeof(Scalar)); 100 /* free up old matrix storage */ 101 FREE(mat->a); if (!mat->singlemalloc) {FREE(mat->i);FREE(mat->j);} 102 aa = mat->a = new_a; ai = mat->i = new_i; aj = mat->j = new_j; 103 mat->singlemalloc = 1; 104 105 rp = aj + ai[row] - 1; ap = aa + ai[row] - 1; 106 rmax = imax[row] = imax[row] + CHUNCKSIZE; 107 mat->mem += CHUNCKSIZE*(sizeof(int) + sizeof(Scalar)); 108 } 109 N = nrow++ - 1; mat->nz++; 110 /* this has too many shifts here; but alternative was slower*/ 111 for ( ii=N; ii>=i; ii-- ) {/* shift values up*/ 112 rp[ii+1] = rp[ii]; 113 ap[ii+1] = ap[ii]; 114 } 115 rp[i] = col; 116 ap[i] = value; 117 noinsert:; 118 a = i + 1; 119 } 120 else { 121 ap[0] = value; rp[0] = col; nrow++; a = 1; 122 } 123 } 124 ailen[row] = nrow; 125 } 126 return 0; 127 } 128 129 static int MatiAIJview(PetscObject obj,Viewer ptr) 130 { 131 Mat aijin = (Mat) obj; 132 Matiaij *aij = (Matiaij *) aijin->data; 133 int i,j; 134 for ( i=0; i<aij->m; i++ ) { 135 printf("row %d:",i); 136 for ( j=aij->i[i]-1; j<aij->i[i+1]-1; j++ ) { 137 #if defined(PETSC_COMPLEX) 138 printf(" %d %g + %g i",aij->j[j]-1,real(aij->a[j]),imag(aij->a[j])); 139 #else 140 printf(" %d %g ",aij->j[j]-1,aij->a[j]); 141 #endif 142 } 143 printf("\n"); 144 } 145 return 0; 146 } 147 148 static int MatiAIJEndAssemble(Mat aijin) 149 { 150 Matiaij *aij = (Matiaij *) aijin->data; 151 int shift = 0,i,j,*ai = aij->i, *aj = aij->j, *imax = aij->imax; 152 int m = aij->m, *ip, N, *ailen = aij->ilen; 153 Scalar *a = aij->a, *ap; 154 155 for ( i=1; i<m; i++ ) { 156 shift += imax[i-1] - ailen[i-1]; 157 if (shift) { 158 ip = aj + ai[i] - 1; ap = a + ai[i] - 1; 159 N = ailen[i]; 160 for ( j=0; j<N; j++ ) { 161 ip[j-shift] = ip[j]; 162 ap[j-shift] = ap[j]; 163 } 164 } 165 ai[i] = ai[i-1] + ailen[i-1]; 166 } 167 shift += imax[m-1] - ailen[m-1]; 168 ai[m] = ai[m-1] + ailen[m-1]; 169 FREE(aij->imax); 170 FREE(aij->ilen); 171 aij->mem -= 2*(aij->m)*sizeof(int); 172 return 0; 173 } 174 175 static int MatiAIJzeroentries(Mat mat) 176 { 177 Matiaij *aij = (Matiaij *) mat->data; 178 Scalar *a = aij->a; 179 int i,n = aij->i[aij->m]-1; 180 for ( i=0; i<n; i++ ) a[i] = 0.0; 181 return 0; 182 183 } 184 static int MatiAIJdestroy(PetscObject obj) 185 { 186 Mat mat = (Mat) obj; 187 Matiaij *aij = (Matiaij *) mat->data; 188 FREE(aij->a); 189 if (!aij->singlemalloc) {FREE(aij->i); FREE(aij->j);} 190 FREE(aij); FREE(mat); 191 return 0; 192 } 193 194 static int MatiAIJCompress(Mat aijin) 195 { 196 return 0; 197 } 198 199 static int MatiAIJinsopt(Mat aijin,int op) 200 { 201 Matiaij *aij = (Matiaij *) aijin->data; 202 if (op == ROW_ORIENTED) aij->roworiented = 1; 203 else if (op == COLUMN_ORIENTED) aij->roworiented = 0; 204 else if (op == COLUMNS_SORTED) aij->sorted = 1; 205 /* doesn't care about sorted rows */ 206 else if (op == NO_NEW_NONZERO_LOCATIONS) aij->nonew = 1; 207 else if (op == YES_NEW_NONZERO_LOCATIONS) aij->nonew = 0; 208 209 if (op == COLUMN_ORIENTED) SETERR(1,"Column oriented input not supported"); 210 return 0; 211 } 212 213 static int MatiAIJgetdiag(Mat aijin,Vec v) 214 { 215 Matiaij *aij = (Matiaij *) aijin->data; 216 int i,j, n; 217 Scalar *x, zero = 0.0; 218 CHKTYPE(v,SEQVECTOR); 219 VecSet(&zero,v); 220 VecGetArray(v,&x); VecGetSize(v,&n); 221 if (n != aij->m) SETERR(1,"Nonconforming matrix and vector"); 222 for ( i=0; i<aij->m; i++ ) { 223 for ( j=aij->i[i]-1; j<aij->i[i+1]-1; j++ ) { 224 if (aij->j[j]-1 == i) { 225 x[i] = aij->a[j]; 226 break; 227 } 228 } 229 } 230 return 0; 231 } 232 233 /* -------------------------------------------------------*/ 234 /* Should check that shapes of vectors and matrices match */ 235 /* -------------------------------------------------------*/ 236 static int MatiAIJmulttrans(Mat aijin,Vec xx,Vec yy) 237 { 238 Matiaij *aij = (Matiaij *) aijin->data; 239 Scalar *x, *y, *v, alpha; 240 int m = aij->m, n, i, *idx; 241 CHKTYPE(xx,SEQVECTOR);CHKTYPE(yy,SEQVECTOR); 242 VecGetArray(xx,&x); VecGetArray(yy,&y); 243 MEMSET(y,0,aij->n*sizeof(Scalar)); 244 y--; /* shift for Fortran start by 1 indexing */ 245 for ( i=0; i<m; i++ ) { 246 idx = aij->j + aij->i[i] - 1; 247 v = aij->a + aij->i[i] - 1; 248 n = aij->i[i+1] - aij->i[i]; 249 alpha = x[i]; 250 /* should inline */ 251 while (n-->0) {y[*idx++] += alpha * *v++;} 252 } 253 return 0; 254 } 255 static int MatiAIJmulttransadd(Mat aijin,Vec xx,Vec zz,Vec yy) 256 { 257 Matiaij *aij = (Matiaij *) aijin->data; 258 Scalar *x, *y, *z, *v, alpha; 259 int m = aij->m, n, i, *idx; 260 CHKTYPE(xx,SEQVECTOR);CHKTYPE(yy,SEQVECTOR); 261 VecGetArray(xx,&x); VecGetArray(yy,&y); 262 if (zz != yy) VecCopy(zz,yy); 263 y--; /* shift for Fortran start by 1 indexing */ 264 for ( i=0; i<m; i++ ) { 265 idx = aij->j + aij->i[i] - 1; 266 v = aij->a + aij->i[i] - 1; 267 n = aij->i[i+1] - aij->i[i]; 268 alpha = x[i]; 269 /* should inline */ 270 while (n-->0) {y[*idx++] += alpha * *v++;} 271 } 272 return 0; 273 } 274 275 static int MatiAIJmult(Mat aijin,Vec xx,Vec yy) 276 { 277 Matiaij *aij = (Matiaij *) aijin->data; 278 Scalar *x, *y, *v, sum; 279 int m = aij->m, n, i, *idx; 280 CHKTYPE(xx,SEQVECTOR);CHKTYPE(yy,SEQVECTOR); 281 VecGetArray(xx,&x); VecGetArray(yy,&y); 282 x--; /* shift for Fortran start by 1 indexing */ 283 for ( i=0; i<m; i++ ) { 284 idx = aij->j + aij->i[i] - 1; 285 v = aij->a + aij->i[i] - 1; 286 n = aij->i[i+1] - aij->i[i]; 287 sum = 0.0; 288 SPARSEDENSEDOT(sum,x,v,idx,n); 289 y[i] = sum; 290 } 291 return 0; 292 } 293 294 static int MatiAIJmultadd(Mat aijin,Vec xx,Vec yy,Vec zz) 295 { 296 Matiaij *aij = (Matiaij *) aijin->data; 297 Scalar *x, *y, *z, *v, sum; 298 int m = aij->m, n, i, *idx; 299 CHKTYPE(xx,SEQVECTOR);CHKTYPE(yy,SEQVECTOR); CHKTYPE(zz,SEQVECTOR); 300 VecGetArray(xx,&x); VecGetArray(yy,&y); VecGetArray(zz,&z); 301 x--; /* shift for Fortran start by 1 indexing */ 302 for ( i=0; i<m; i++ ) { 303 idx = aij->j + aij->i[i] - 1; 304 v = aij->a + aij->i[i] - 1; 305 n = aij->i[i+1] - aij->i[i]; 306 sum = y[i]; 307 SPARSEDENSEDOT(sum,x,v,idx,n); 308 y[i] = sum; 309 } 310 return 0; 311 } 312 313 /* 314 Adds diagonal pointers to sparse matrix structure. 315 */ 316 317 static int MatiAIJmarkdiag(Matiaij *aij) 318 { 319 int i,j, n, *diag;; 320 diag = (int *) MALLOC( aij->m*sizeof(int)); CHKPTR(diag); 321 for ( i=0; i<aij->m; i++ ) { 322 for ( j=aij->i[i]-1; j<aij->i[i+1]-1; j++ ) { 323 if (aij->j[j]-1 == i) { 324 diag[i] = j + 1; 325 break; 326 } 327 } 328 } 329 aij->diag = diag; 330 return 0; 331 } 332 333 static int MatiAIJrelax(Mat matin,Vec bb,double omega,int flag,IS is, 334 int its,Vec xx) 335 { 336 Matiaij *mat = (Matiaij *) matin->data; 337 Scalar *x, *b, zero = 0.0, d, *xs, sum, *v = mat->a; 338 int ierr,one = 1, tmp, *idx, *diag; 339 int n = mat->n, m = mat->m, i, j; 340 341 if (is) SETERR(1,"No support for ordering in relaxation"); 342 if (flag & SOR_ZERO_INITIAL_GUESS) { 343 if (ierr = VecSet(&zero,xx)) return ierr; 344 } 345 VecGetArray(xx,&x); VecGetArray(bb,&b); 346 if (!mat->diag) {if (ierr = MatiAIJmarkdiag(mat)) return ierr;} 347 diag = mat->diag; 348 xs = x - 1; /* shifted by one for index start of a or mat->j*/ 349 while (its--) { 350 if (flag & SOR_FORWARD_SWEEP){ 351 for ( i=0; i<m; i++ ) { 352 d = mat->a[diag[i]-1]; 353 n = mat->i[i+1] - mat->i[i]; 354 idx = mat->j + mat->i[i] - 1; 355 v = mat->a + mat->i[i] - 1; 356 sum = b[i]; 357 SPARSEDENSEMDOT(sum,xs,v,idx,n); 358 x[i] = (1. - omega)*x[i] + omega*(sum/d + x[i]); 359 } 360 } 361 if (flag & SOR_BACKWARD_SWEEP){ 362 for ( i=m-1; i>=0; i-- ) { 363 d = mat->a[diag[i] - 1]; 364 n = mat->i[i+1] - mat->i[i]; 365 idx = mat->j + mat->i[i] - 1; 366 v = mat->a + mat->i[i] - 1; 367 sum = b[i]; 368 SPARSEDENSEMDOT(sum,xs,v,idx,n); 369 x[i] = (1. - omega)*x[i] + omega*(sum/d + x[i]); 370 } 371 } 372 } 373 return 0; 374 } 375 376 static int MatiAIJnz(Mat matin,int *nz) 377 { 378 Matiaij *mat = (Matiaij *) matin->data; 379 *nz = mat->nz; 380 return 0; 381 } 382 static int MatiAIJmem(Mat matin,int *nz) 383 { 384 Matiaij *mat = (Matiaij *) matin->data; 385 *nz = mat->mem; 386 return 0; 387 } 388 389 int MatiAIJLUFactorSymbolic(Mat,IS,IS,Mat*); 390 int MatiAIJLUFactorNumeric(Mat,Mat); 391 int MatiAIJSolve(Mat,Vec,Vec); 392 393 /* -------------------------------------------------------------------*/ 394 static struct _MatOps MatOps = {MatiAIJAddValues,MatiAIJAddValues, 395 0, 0, 396 MatiAIJmult,MatiAIJmultadd,MatiAIJmulttrans,MatiAIJmulttransadd, 397 MatiAIJSolve,0,0,0, 398 0,0, 399 MatiAIJrelax, 400 0, 401 MatiAIJnz,MatiAIJmem,0, 402 0, 403 MatiAIJgetdiag,0,0, 404 0,MatiAIJEndAssemble, 405 MatiAIJCompress, 406 MatiAIJinsopt,MatiAIJzeroentries,0,MatAIJreorder, 407 MatiAIJLUFactorSymbolic,MatiAIJLUFactorNumeric,0,0 }; 408 409 410 411 /*@ 412 413 MatCreateSequentialAIJ - Creates a sparse matrix in AIJ format. 414 415 Input Parameters: 416 . m,n - number of rows and columns 417 . nz - total number nonzeros in matrix 418 . nzz - number of nonzeros per row or null 419 . You must leave room for the diagonal entry even if it is zero. 420 421 Output parameters: 422 . newmat - the matrix 423 424 Keywords: matrix, aij, compressed row, sparse 425 @*/ 426 int MatCreateSequentialAIJ(int m,int n,int nz,int *nnz, Mat *newmat) 427 { 428 Mat mat; 429 Matiaij *aij; 430 int i,rl,len; 431 *newmat = 0; 432 CREATEHEADER(mat,_Mat); 433 mat->data = (void *) (aij = NEW(Matiaij)); CHKPTR(aij); 434 mat->cookie = MAT_COOKIE; 435 mat->ops = &MatOps; 436 mat->destroy = MatiAIJdestroy; 437 mat->view = MatiAIJview; 438 mat->type = MATAIJSEQ; 439 mat->factor = 0; 440 mat->row = 0; 441 mat->col = 0; 442 aij->m = m; 443 aij->n = n; 444 aij->imax = (int *) MALLOC( m*sizeof(int) ); CHKPTR(aij->imax); 445 aij->mem = m*sizeof(int) + sizeof(Matiaij); 446 if (!nnz) { 447 if (nz <= 0) nz = 1; 448 for ( i=0; i<m; i++ ) aij->imax[i] = nz; 449 nz = nz*m; 450 } 451 else { 452 nz = 0; 453 for ( i=0; i<m; i++ ) {aij->imax[i] = nnz[i]; nz += nnz[i];} 454 } 455 456 /* allocate the matrix space */ 457 aij->nz = nz; 458 len = nz*(sizeof(int) + sizeof(Scalar)) + (aij->m+1)*sizeof(int); 459 aij->a = (Scalar *) MALLOC( len ); CHKPTR(aij->a); 460 aij->j = (int *) (aij->a + nz); 461 aij->i = aij->j + nz; 462 aij->singlemalloc = 1; 463 aij->mem += len; 464 465 aij->i[0] = 1; 466 for (i=1; i<m+1; i++) { 467 aij->i[i] = aij->i[i-1] + aij->imax[i-1]; 468 } 469 470 /* aij->ilen will count nonzeros in each row so far. */ 471 aij->ilen = (int *) MALLOC((aij->m)*sizeof(int)); 472 473 /* stick in zeros along diagonal */ 474 for ( i=0; i<aij->m; i++ ) { 475 aij->ilen[i] = 1; 476 aij->j[aij->i[i]-1] = i+1; 477 aij->a[aij->i[i]-1] = 0.0; 478 } 479 aij->nz = aij->m; 480 aij->mem += (aij->m)*sizeof(int) + len; 481 482 aij->sorted = 0; 483 aij->roworiented = 1; 484 aij->nonew = 0; 485 aij->diag = 0; 486 487 *newmat = mat; 488 return 0; 489 } 490