1 2 #ifdef PETSC_RCS_HEADER 3 static char vcid[] = "$Id: aij.c,v 1.234 1997/08/28 14:23:17 curfman Exp bsmith $"; 4 #endif 5 6 /* 7 Defines the basic matrix operations for the AIJ (compressed row) 8 matrix storage format. 9 */ 10 11 #include "pinclude/pviewer.h" 12 #include "sys.h" 13 #include "src/mat/impls/aij/seq/aij.h" 14 #include "src/vec/vecimpl.h" 15 #include "src/inline/spops.h" 16 #include "src/inline/dot.h" 17 #include "src/inline/bitarray.h" 18 19 /* 20 Basic AIJ format ILU based on drop tolerance 21 */ 22 #undef __FUNC__ 23 #define __FUNC__ "MatILUDTFactor_SeqAIJ" 24 int MatILUDTFactor_SeqAIJ(Mat A,double dt,int maxnz,IS row,IS col,Mat *fact) 25 { 26 /* Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; */ 27 int ierr = 1; 28 29 SETERRQ(ierr,0,"Not implemented"); 30 } 31 32 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 33 34 #undef __FUNC__ 35 #define __FUNC__ "MatGetRowIJ_SeqAIJ" 36 int MatGetRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *m,int **ia,int **ja, 37 PetscTruth *done) 38 { 39 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 40 int ierr,i,ishift; 41 42 *m = A->m; 43 if (!ia) return 0; 44 ishift = a->indexshift; 45 if (symmetric) { 46 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja); CHKERRQ(ierr); 47 } else if (oshift == 0 && ishift == -1) { 48 int nz = a->i[a->m]; 49 /* malloc space and subtract 1 from i and j indices */ 50 *ia = (int *) PetscMalloc( (a->m+1)*sizeof(int) ); CHKPTRQ(*ia); 51 *ja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(*ja); 52 for ( i=0; i<nz; i++ ) (*ja)[i] = a->j[i] - 1; 53 for ( i=0; i<a->m+1; i++ ) (*ia)[i] = a->i[i] - 1; 54 } else if (oshift == 1 && ishift == 0) { 55 int nz = a->i[a->m] + 1; 56 /* malloc space and add 1 to i and j indices */ 57 *ia = (int *) PetscMalloc( (a->m+1)*sizeof(int) ); CHKPTRQ(*ia); 58 *ja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(*ja); 59 for ( i=0; i<nz; i++ ) (*ja)[i] = a->j[i] + 1; 60 for ( i=0; i<a->m+1; i++ ) (*ia)[i] = a->i[i] + 1; 61 } else { 62 *ia = a->i; *ja = a->j; 63 } 64 65 return 0; 66 } 67 68 #undef __FUNC__ 69 #define __FUNC__ "MatRestoreRowIJ_SeqAIJ" 70 int MatRestoreRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia,int **ja, 71 PetscTruth *done) 72 { 73 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 74 int ishift = a->indexshift; 75 76 if (!ia) return 0; 77 if (symmetric || (oshift == 0 && ishift == -1) || (oshift == 1 && ishift == 0)) { 78 PetscFree(*ia); 79 PetscFree(*ja); 80 } 81 return 0; 82 } 83 84 #undef __FUNC__ 85 #define __FUNC__ "MatGetColumnIJ_SeqAIJ" 86 int MatGetColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 87 PetscTruth *done) 88 { 89 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 90 int ierr,i,ishift = a->indexshift,*collengths,*cia,*cja,n = A->n,m = A->m; 91 int nz = a->i[m]+ishift,row,*jj,mr,col; 92 93 *nn = A->n; 94 if (!ia) return 0; 95 if (symmetric) { 96 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja); CHKERRQ(ierr); 97 } else { 98 collengths = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(collengths); 99 PetscMemzero(collengths,n*sizeof(int)); 100 cia = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(cia); 101 cja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(cja); 102 jj = a->j; 103 for ( i=0; i<nz; i++ ) { 104 collengths[jj[i] + ishift]++; 105 } 106 cia[0] = oshift; 107 for ( i=0; i<n; i++) { 108 cia[i+1] = cia[i] + collengths[i]; 109 } 110 PetscMemzero(collengths,n*sizeof(int)); 111 jj = a->j; 112 for ( row=0; row<m; row++ ) { 113 mr = a->i[row+1] - a->i[row]; 114 for ( i=0; i<mr; i++ ) { 115 col = *jj++ + ishift; 116 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 117 } 118 } 119 PetscFree(collengths); 120 *ia = cia; *ja = cja; 121 } 122 123 return 0; 124 } 125 126 #undef __FUNC__ 127 #define __FUNC__ "MatRestoreColumnIJ_SeqAIJ" 128 int MatRestoreColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia, 129 int **ja,PetscTruth *done) 130 { 131 if (!ia) return 0; 132 133 PetscFree(*ia); 134 PetscFree(*ja); 135 136 return 0; 137 } 138 139 #define CHUNKSIZE 15 140 141 #undef __FUNC__ 142 #define __FUNC__ "MatSetValues_SeqAIJ" 143 int MatSetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 144 { 145 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 146 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax, N, sorted = a->sorted; 147 int *imax = a->imax, *ai = a->i, *ailen = a->ilen,roworiented = a->roworiented; 148 int *aj = a->j, nonew = a->nonew,shift = a->indexshift; 149 Scalar *ap,value, *aa = a->a; 150 151 for ( k=0; k<m; k++ ) { /* loop over added rows */ 152 row = im[k]; 153 #if defined(PETSC_BOPT_g) 154 if (row < 0) SETERRQ(1,0,"Negative row"); 155 if (row >= a->m) SETERRQ(1,0,"Row too large"); 156 #endif 157 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 158 rmax = imax[row]; nrow = ailen[row]; 159 low = 0; 160 for ( l=0; l<n; l++ ) { /* loop over added columns */ 161 #if defined(PETSC_BOPT_g) 162 if (in[l] < 0) SETERRQ(1,0,"Negative column"); 163 if (in[l] >= a->n) SETERRQ(1,0,"Column too large"); 164 #endif 165 col = in[l] - shift; 166 if (roworiented) { 167 value = *v++; 168 } 169 else { 170 value = v[k + l*m]; 171 } 172 if (!sorted) low = 0; high = nrow; 173 while (high-low > 5) { 174 t = (low+high)/2; 175 if (rp[t] > col) high = t; 176 else low = t; 177 } 178 for ( i=low; i<high; i++ ) { 179 if (rp[i] > col) break; 180 if (rp[i] == col) { 181 if (is == ADD_VALUES) ap[i] += value; 182 else ap[i] = value; 183 goto noinsert; 184 } 185 } 186 if (nonew == 1) goto noinsert; 187 else if (nonew == -1) SETERRQ(1,0,"Inserting a new nonzero in the matrix"); 188 if (nrow >= rmax) { 189 /* there is no extra room in row, therefore enlarge */ 190 int new_nz = ai[a->m] + CHUNKSIZE,len,*new_i,*new_j; 191 Scalar *new_a; 192 193 if (nonew == -2) SETERRQ(1,0,"Inserting a new nonzero in the matrix"); 194 195 /* malloc new storage space */ 196 len = new_nz*(sizeof(int)+sizeof(Scalar))+(a->m+1)*sizeof(int); 197 new_a = (Scalar *) PetscMalloc( len ); CHKPTRQ(new_a); 198 new_j = (int *) (new_a + new_nz); 199 new_i = new_j + new_nz; 200 201 /* copy over old data into new slots */ 202 for ( ii=0; ii<row+1; ii++ ) {new_i[ii] = ai[ii];} 203 for ( ii=row+1; ii<a->m+1; ii++ ) {new_i[ii] = ai[ii]+CHUNKSIZE;} 204 PetscMemcpy(new_j,aj,(ai[row]+nrow+shift)*sizeof(int)); 205 len = (new_nz - CHUNKSIZE - ai[row] - nrow - shift); 206 PetscMemcpy(new_j+ai[row]+shift+nrow+CHUNKSIZE,aj+ai[row]+shift+nrow, 207 len*sizeof(int)); 208 PetscMemcpy(new_a,aa,(ai[row]+nrow+shift)*sizeof(Scalar)); 209 PetscMemcpy(new_a+ai[row]+shift+nrow+CHUNKSIZE,aa+ai[row]+shift+nrow, 210 len*sizeof(Scalar)); 211 /* free up old matrix storage */ 212 PetscFree(a->a); 213 if (!a->singlemalloc) {PetscFree(a->i);PetscFree(a->j);} 214 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 215 a->singlemalloc = 1; 216 217 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 218 rmax = imax[row] = imax[row] + CHUNKSIZE; 219 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + sizeof(Scalar))); 220 a->maxnz += CHUNKSIZE; 221 a->reallocs++; 222 } 223 N = nrow++ - 1; a->nz++; 224 /* shift up all the later entries in this row */ 225 for ( ii=N; ii>=i; ii-- ) { 226 rp[ii+1] = rp[ii]; 227 ap[ii+1] = ap[ii]; 228 } 229 rp[i] = col; 230 ap[i] = value; 231 noinsert:; 232 low = i + 1; 233 } 234 ailen[row] = nrow; 235 } 236 return 0; 237 } 238 239 #undef __FUNC__ 240 #define __FUNC__ "MatGetValues_SeqAIJ" 241 int MatGetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 242 { 243 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 244 int *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j; 245 int *ai = a->i, *ailen = a->ilen, shift = a->indexshift; 246 Scalar *ap, *aa = a->a, zero = 0.0; 247 248 for ( k=0; k<m; k++ ) { /* loop over rows */ 249 row = im[k]; 250 if (row < 0) SETERRQ(1,0,"Negative row"); 251 if (row >= a->m) SETERRQ(1,0,"Row too large"); 252 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 253 nrow = ailen[row]; 254 for ( l=0; l<n; l++ ) { /* loop over columns */ 255 if (in[l] < 0) SETERRQ(1,0,"Negative column"); 256 if (in[l] >= a->n) SETERRQ(1,0,"Column too large"); 257 col = in[l] - shift; 258 high = nrow; low = 0; /* assume unsorted */ 259 while (high-low > 5) { 260 t = (low+high)/2; 261 if (rp[t] > col) high = t; 262 else low = t; 263 } 264 for ( i=low; i<high; i++ ) { 265 if (rp[i] > col) break; 266 if (rp[i] == col) { 267 *v++ = ap[i]; 268 goto finished; 269 } 270 } 271 *v++ = zero; 272 finished:; 273 } 274 } 275 return 0; 276 } 277 278 279 #undef __FUNC__ 280 #define __FUNC__ "MatView_SeqAIJ_Binary" 281 extern int MatView_SeqAIJ_Binary(Mat A,Viewer viewer) 282 { 283 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 284 int i, fd, *col_lens, ierr; 285 286 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 287 col_lens = (int *) PetscMalloc( (4+a->m)*sizeof(int) ); CHKPTRQ(col_lens); 288 col_lens[0] = MAT_COOKIE; 289 col_lens[1] = a->m; 290 col_lens[2] = a->n; 291 col_lens[3] = a->nz; 292 293 /* store lengths of each row and write (including header) to file */ 294 for ( i=0; i<a->m; i++ ) { 295 col_lens[4+i] = a->i[i+1] - a->i[i]; 296 } 297 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,BINARY_INT,1); CHKERRQ(ierr); 298 PetscFree(col_lens); 299 300 /* store column indices (zero start index) */ 301 if (a->indexshift) { 302 for ( i=0; i<a->nz; i++ ) a->j[i]--; 303 } 304 ierr = PetscBinaryWrite(fd,a->j,a->nz,BINARY_INT,0); CHKERRQ(ierr); 305 if (a->indexshift) { 306 for ( i=0; i<a->nz; i++ ) a->j[i]++; 307 } 308 309 /* store nonzero values */ 310 ierr = PetscBinaryWrite(fd,a->a,a->nz,BINARY_SCALAR,0); CHKERRQ(ierr); 311 return 0; 312 } 313 314 #undef __FUNC__ 315 #define __FUNC__ "MatView_SeqAIJ_ASCII" 316 extern int MatView_SeqAIJ_ASCII(Mat A,Viewer viewer) 317 { 318 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 319 int ierr, i,j, m = a->m, shift = a->indexshift, format, flg1,flg2; 320 FILE *fd; 321 char *outputname; 322 323 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 324 ierr = ViewerFileGetOutputname_Private(viewer,&outputname); CHKERRQ(ierr); 325 ierr = ViewerGetFormat(viewer,&format); 326 if (format == VIEWER_FORMAT_ASCII_INFO) { 327 return 0; 328 } 329 else if (format == VIEWER_FORMAT_ASCII_INFO_LONG) { 330 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_no_inode",&flg1); CHKERRQ(ierr); 331 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg2); CHKERRQ(ierr); 332 if (flg1 || flg2) fprintf(fd," not using I-node routines\n"); 333 else fprintf(fd," using I-node routines: found %d nodes, limit used is %d\n", 334 a->inode.node_count,a->inode.limit); 335 } 336 else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 337 int nofinalvalue = 0; 338 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != a->n-!shift)) { 339 nofinalvalue = 1; 340 } 341 fprintf(fd,"%% Size = %d %d \n",m,a->n); 342 fprintf(fd,"%% Nonzeros = %d \n",a->nz); 343 fprintf(fd,"zzz = zeros(%d,3);\n",a->nz+nofinalvalue); 344 fprintf(fd,"zzz = [\n"); 345 346 for (i=0; i<m; i++) { 347 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 348 #if defined(PETSC_COMPLEX) 349 fprintf(fd,"%d %d %18.16e + %18.16e i \n",i+1,a->j[j]+!shift,real(a->a[j]), 350 imag(a->a[j])); 351 #else 352 fprintf(fd,"%d %d %18.16e\n", i+1, a->j[j]+!shift, a->a[j]); 353 #endif 354 } 355 } 356 if (nofinalvalue) { 357 fprintf(fd,"%d %d %18.16e\n", m, a->n, 0.0); 358 } 359 fprintf(fd,"];\n %s = spconvert(zzz);\n",outputname); 360 } 361 else if (format == VIEWER_FORMAT_ASCII_COMMON) { 362 for ( i=0; i<m; i++ ) { 363 fprintf(fd,"row %d:",i); 364 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 365 #if defined(PETSC_COMPLEX) 366 if (imag(a->a[j]) > 0.0 && real(a->a[j]) != 0.0) 367 fprintf(fd," %d %g + %g i",a->j[j]+shift,real(a->a[j]),imag(a->a[j])); 368 else if (imag(a->a[j]) < 0.0 && real(a->a[j]) != 0.0) 369 fprintf(fd," %d %g - %g i",a->j[j]+shift,real(a->a[j]),-imag(a->a[j])); 370 else if (real(a->a[j]) != 0.0) 371 fprintf(fd," %d %g ",a->j[j]+shift,real(a->a[j])); 372 #else 373 if (a->a[j] != 0.0) fprintf(fd," %d %g ",a->j[j]+shift,a->a[j]); 374 #endif 375 } 376 fprintf(fd,"\n"); 377 } 378 } 379 else if (format == VIEWER_FORMAT_ASCII_SYMMODU) { 380 int nzd=0, fshift=1, *sptr; 381 sptr = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(sptr); 382 for ( i=0; i<m; i++ ) { 383 sptr[i] = nzd+1; 384 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 385 if (a->j[j] >= i) { 386 #if defined(PETSC_COMPLEX) 387 if (imag(a->a[j]) != 0.0 || real(a->a[j]) != 0.0) nzd++; 388 #else 389 if (a->a[j] != 0.0) nzd++; 390 #endif 391 } 392 } 393 } 394 sptr[m] = nzd+1; 395 fprintf(fd," %d %d\n\n",m,nzd); 396 for ( i=0; i<m+1; i+=6 ) { 397 if (i+4<m) fprintf(fd," %d %d %d %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5]); 398 else if (i+3<m) fprintf(fd," %d %d %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4]); 399 else if (i+2<m) fprintf(fd," %d %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]); 400 else if (i+1<m) fprintf(fd," %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2]); 401 else if (i<m) fprintf(fd," %d %d\n",sptr[i],sptr[i+1]); 402 else fprintf(fd," %d\n",sptr[i]); 403 } 404 fprintf(fd,"\n"); 405 PetscFree(sptr); 406 for ( i=0; i<m; i++ ) { 407 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 408 if (a->j[j] >= i) fprintf(fd," %d ",a->j[j]+fshift); 409 } 410 fprintf(fd,"\n"); 411 } 412 fprintf(fd,"\n"); 413 for ( i=0; i<m; i++ ) { 414 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 415 if (a->j[j] >= i) { 416 #if defined(PETSC_COMPLEX) 417 if (imag(a->a[j]) != 0.0 || real(a->a[j]) != 0.0) 418 fprintf(fd," %18.16e %18.16e ",real(a->a[j]),imag(a->a[j])); 419 #else 420 if (a->a[j] != 0.0) fprintf(fd," %18.16e ",a->a[j]); 421 #endif 422 } 423 } 424 fprintf(fd,"\n"); 425 } 426 } 427 else { 428 for ( i=0; i<m; i++ ) { 429 fprintf(fd,"row %d:",i); 430 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 431 #if defined(PETSC_COMPLEX) 432 if (imag(a->a[j]) > 0.0) { 433 fprintf(fd," %d %g + %g i",a->j[j]+shift,real(a->a[j]),imag(a->a[j])); 434 } else if (imag(a->a[j]) < 0.0) { 435 fprintf(fd," %d %g - %g i",a->j[j]+shift,real(a->a[j]),-imag(a->a[j])); 436 } 437 else { 438 fprintf(fd," %d %g ",a->j[j]+shift,real(a->a[j])); 439 } 440 #else 441 fprintf(fd," %d %g ",a->j[j]+shift,a->a[j]); 442 #endif 443 } 444 fprintf(fd,"\n"); 445 } 446 } 447 fflush(fd); 448 return 0; 449 } 450 451 #undef __FUNC__ 452 #define __FUNC__ "MatView_SeqAIJ_Draw" 453 extern int MatView_SeqAIJ_Draw(Mat A,Viewer viewer) 454 { 455 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 456 int ierr, i,j, m = a->m, shift = a->indexshift,pause,color; 457 int format; 458 double xl,yl,xr,yr,w,h,xc,yc,scale = 1.0,x_l,x_r,y_l,y_r,maxv = 0.0; 459 Draw draw; 460 DrawButton button; 461 PetscTruth isnull; 462 463 ierr = ViewerDrawGetDraw(viewer,&draw); CHKERRQ(ierr); 464 ierr = DrawClear(draw); CHKERRQ(ierr); 465 ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr); 466 ierr = DrawIsNull(draw,&isnull); CHKERRQ(ierr); if (isnull) return 0; 467 468 xr = a->n; yr = a->m; h = yr/10.0; w = xr/10.0; 469 xr += w; yr += h; xl = -w; yl = -h; 470 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 471 /* loop over matrix elements drawing boxes */ 472 473 if (format != VIEWER_FORMAT_DRAW_CONTOUR) { 474 /* Blue for negative, Cyan for zero and Red for positive */ 475 color = DRAW_BLUE; 476 for ( i=0; i<m; i++ ) { 477 y_l = m - i - 1.0; y_r = y_l + 1.0; 478 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 479 x_l = a->j[j] + shift; x_r = x_l + 1.0; 480 #if defined(PETSC_COMPLEX) 481 if (real(a->a[j]) >= 0.) continue; 482 #else 483 if (a->a[j] >= 0.) continue; 484 #endif 485 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 486 } 487 } 488 color = DRAW_CYAN; 489 for ( i=0; i<m; i++ ) { 490 y_l = m - i - 1.0; y_r = y_l + 1.0; 491 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 492 x_l = a->j[j] + shift; x_r = x_l + 1.0; 493 if (a->a[j] != 0.) continue; 494 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 495 } 496 } 497 color = DRAW_RED; 498 for ( i=0; i<m; i++ ) { 499 y_l = m - i - 1.0; y_r = y_l + 1.0; 500 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 501 x_l = a->j[j] + shift; x_r = x_l + 1.0; 502 #if defined(PETSC_COMPLEX) 503 if (real(a->a[j]) <= 0.) continue; 504 #else 505 if (a->a[j] <= 0.) continue; 506 #endif 507 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 508 } 509 } 510 } else { 511 /* use contour shading to indicate magnitude of values */ 512 /* first determine max of all nonzero values */ 513 int nz = a->nz,count; 514 Draw popup; 515 516 for ( i=0; i<nz; i++ ) { 517 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 518 } 519 ierr = DrawCreatePopUp(draw,&popup); CHKERRQ(ierr); 520 ierr = DrawScalePopup(popup,0.0,maxv); CHKERRQ(ierr); 521 count = 0; 522 for ( i=0; i<m; i++ ) { 523 y_l = m - i - 1.0; y_r = y_l + 1.0; 524 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 525 x_l = a->j[j] + shift; x_r = x_l + 1.0; 526 color = 32 + (int) ((200.0 - 32.0)*PetscAbsScalar(a->a[count])/maxv); 527 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 528 count++; 529 } 530 } 531 } 532 DrawFlush(draw); 533 DrawGetPause(draw,&pause); 534 if (pause >= 0) { PetscSleep(pause); return 0;} 535 536 /* allow the matrix to zoom or shrink */ 537 ierr = DrawCheckResizedWindow(draw); 538 ierr = DrawGetMouseButton(draw,&button,&xc,&yc,0,0); 539 while (button != BUTTON_RIGHT) { 540 DrawClear(draw); 541 if (button == BUTTON_LEFT) scale = .5; 542 else if (button == BUTTON_CENTER) scale = 2.; 543 xl = scale*(xl + w - xc) + xc - w*scale; 544 xr = scale*(xr - w - xc) + xc + w*scale; 545 yl = scale*(yl + h - yc) + yc - h*scale; 546 yr = scale*(yr - h - yc) + yc + h*scale; 547 w *= scale; h *= scale; 548 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 549 if (format != VIEWER_FORMAT_DRAW_CONTOUR) { 550 /* Blue for negative, Cyan for zero and Red for positive */ 551 color = DRAW_BLUE; 552 for ( i=0; i<m; i++ ) { 553 y_l = m - i - 1.0; y_r = y_l + 1.0; 554 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 555 x_l = a->j[j] + shift; x_r = x_l + 1.0; 556 #if defined(PETSC_COMPLEX) 557 if (real(a->a[j]) >= 0.) continue; 558 #else 559 if (a->a[j] >= 0.) continue; 560 #endif 561 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 562 } 563 } 564 color = DRAW_CYAN; 565 for ( i=0; i<m; i++ ) { 566 y_l = m - i - 1.0; y_r = y_l + 1.0; 567 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 568 x_l = a->j[j] + shift; x_r = x_l + 1.0; 569 if (a->a[j] != 0.) continue; 570 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 571 } 572 } 573 color = DRAW_RED; 574 for ( i=0; i<m; i++ ) { 575 y_l = m - i - 1.0; y_r = y_l + 1.0; 576 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 577 x_l = a->j[j] + shift; x_r = x_l + 1.0; 578 #if defined(PETSC_COMPLEX) 579 if (real(a->a[j]) <= 0.) continue; 580 #else 581 if (a->a[j] <= 0.) continue; 582 #endif 583 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 584 } 585 } 586 } else { 587 /* use contour shading to indicate magnitude of values */ 588 int count = 0; 589 for ( i=0; i<m; i++ ) { 590 y_l = m - i - 1.0; y_r = y_l + 1.0; 591 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 592 x_l = a->j[j] + shift; x_r = x_l + 1.0; 593 color = 32 + (int) ((200.0 - 32.0)*PetscAbsScalar(a->a[count])/maxv); 594 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 595 count++; 596 } 597 } 598 } 599 600 ierr = DrawCheckResizedWindow(draw); 601 ierr = DrawGetMouseButton(draw,&button,&xc,&yc,0,0); 602 } 603 return 0; 604 } 605 606 #undef __FUNC__ 607 #define __FUNC__ "MatView_SeqAIJ" 608 int MatView_SeqAIJ(PetscObject obj,Viewer viewer) 609 { 610 Mat A = (Mat) obj; 611 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 612 ViewerType vtype; 613 int ierr; 614 615 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 616 if (vtype == MATLAB_VIEWER) { 617 return ViewerMatlabPutSparse_Private(viewer,a->m,a->n,a->nz,a->a,a->i,a->j); 618 } 619 else if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER){ 620 return MatView_SeqAIJ_ASCII(A,viewer); 621 } 622 else if (vtype == BINARY_FILE_VIEWER) { 623 return MatView_SeqAIJ_Binary(A,viewer); 624 } 625 else if (vtype == DRAW_VIEWER) { 626 return MatView_SeqAIJ_Draw(A,viewer); 627 } 628 return 0; 629 } 630 631 extern int Mat_AIJ_CheckInode(Mat); 632 #undef __FUNC__ 633 #define __FUNC__ "MatAssemblyEnd_SeqAIJ" 634 int MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 635 { 636 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 637 int fshift = 0,i,j,*ai = a->i, *aj = a->j, *imax = a->imax,ierr; 638 int m = a->m, *ip, N, *ailen = a->ilen,shift = a->indexshift,rmax = 0; 639 Scalar *aa = a->a, *ap; 640 641 if (mode == MAT_FLUSH_ASSEMBLY) return 0; 642 643 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 644 for ( i=1; i<m; i++ ) { 645 /* move each row back by the amount of empty slots (fshift) before it*/ 646 fshift += imax[i-1] - ailen[i-1]; 647 rmax = PetscMax(rmax,ailen[i]); 648 if (fshift) { 649 ip = aj + ai[i] + shift; ap = aa + ai[i] + shift; 650 N = ailen[i]; 651 for ( j=0; j<N; j++ ) { 652 ip[j-fshift] = ip[j]; 653 ap[j-fshift] = ap[j]; 654 } 655 } 656 ai[i] = ai[i-1] + ailen[i-1]; 657 } 658 if (m) { 659 fshift += imax[m-1] - ailen[m-1]; 660 ai[m] = ai[m-1] + ailen[m-1]; 661 } 662 /* reset ilen and imax for each row */ 663 for ( i=0; i<m; i++ ) { 664 ailen[i] = imax[i] = ai[i+1] - ai[i]; 665 } 666 a->nz = ai[m] + shift; 667 668 /* diagonals may have moved, so kill the diagonal pointers */ 669 if (fshift && a->diag) { 670 PetscFree(a->diag); 671 PLogObjectMemory(A,-(m+1)*sizeof(int)); 672 a->diag = 0; 673 } 674 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %d X %d; storage space: %d unneeded, %d used\n", 675 m,a->n,fshift,a->nz); 676 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues is %d\n", 677 a->reallocs); 678 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Most nonzeros in any row is %d\n",rmax); 679 a->reallocs = 0; 680 A->info.nz_unneeded = (double)fshift; 681 682 /* check out for identical nodes. If found, use inode functions */ 683 ierr = Mat_AIJ_CheckInode(A); CHKERRQ(ierr); 684 return 0; 685 } 686 687 #undef __FUNC__ 688 #define __FUNC__ "MatZeroEntries_SeqAIJ" 689 int MatZeroEntries_SeqAIJ(Mat A) 690 { 691 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 692 PetscMemzero(a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar)); 693 return 0; 694 } 695 696 #undef __FUNC__ 697 #define __FUNC__ "MatDestroy_SeqAIJ" 698 int MatDestroy_SeqAIJ(PetscObject obj) 699 { 700 Mat A = (Mat) obj; 701 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 702 703 #if defined(PETSC_LOG) 704 PLogObjectState(obj,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 705 #endif 706 PetscFree(a->a); 707 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 708 if (a->diag) PetscFree(a->diag); 709 if (a->ilen) PetscFree(a->ilen); 710 if (a->imax) PetscFree(a->imax); 711 if (a->solve_work) PetscFree(a->solve_work); 712 if (a->inode.size) PetscFree(a->inode.size); 713 PetscFree(a); 714 715 PLogObjectDestroy(A); 716 PetscHeaderDestroy(A); 717 return 0; 718 } 719 720 #undef __FUNC__ 721 #define __FUNC__ "MatCompress_SeqAIJ" 722 int MatCompress_SeqAIJ(Mat A) 723 { 724 return 0; 725 } 726 727 #undef __FUNC__ 728 #define __FUNC__ "MatSetOption_SeqAIJ" 729 int MatSetOption_SeqAIJ(Mat A,MatOption op) 730 { 731 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 732 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 733 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 734 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 735 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 736 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 737 else if (op == MAT_NEW_NONZERO_LOCATION_ERROR) a->nonew = -1; 738 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERROR) a->nonew = -2; 739 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 740 else if (op == MAT_ROWS_SORTED || 741 op == MAT_ROWS_UNSORTED || 742 op == MAT_SYMMETRIC || 743 op == MAT_STRUCTURALLY_SYMMETRIC || 744 op == MAT_YES_NEW_DIAGONALS || 745 op == MAT_IGNORE_OFF_PROC_ENTRIES) 746 PLogInfo(A,"Info:MatSetOption_SeqAIJ:Option ignored\n"); 747 else if (op == MAT_NO_NEW_DIAGONALS) 748 {SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS");} 749 else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 750 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 751 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 752 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 753 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 754 else 755 {SETERRQ(PETSC_ERR_SUP,0,"unknown option");} 756 return 0; 757 } 758 759 #undef __FUNC__ 760 #define __FUNC__ "MatGetDiagonal_SeqAIJ" 761 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 762 { 763 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 764 int i,j, n,shift = a->indexshift; 765 Scalar *x, zero = 0.0; 766 767 VecSet(&zero,v); 768 VecGetArray_Fast(v,x); VecGetLocalSize(v,&n); 769 if (n != a->m) SETERRQ(1,0,"Nonconforming matrix and vector"); 770 for ( i=0; i<a->m; i++ ) { 771 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 772 if (a->j[j]+shift == i) { 773 x[i] = a->a[j]; 774 break; 775 } 776 } 777 } 778 return 0; 779 } 780 781 /* -------------------------------------------------------*/ 782 /* Should check that shapes of vectors and matrices match */ 783 /* -------------------------------------------------------*/ 784 #undef __FUNC__ 785 #define __FUNC__ "MatMultTrans_SeqAIJ" 786 int MatMultTrans_SeqAIJ(Mat A,Vec xx,Vec yy) 787 { 788 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 789 Scalar *x, *y, *v, alpha; 790 int m = a->m, n, i, *idx, shift = a->indexshift; 791 792 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 793 PetscMemzero(y,a->n*sizeof(Scalar)); 794 y = y + shift; /* shift for Fortran start by 1 indexing */ 795 for ( i=0; i<m; i++ ) { 796 idx = a->j + a->i[i] + shift; 797 v = a->a + a->i[i] + shift; 798 n = a->i[i+1] - a->i[i]; 799 alpha = x[i]; 800 while (n-->0) {y[*idx++] += alpha * *v++;} 801 } 802 PLogFlops(2*a->nz - a->n); 803 return 0; 804 } 805 806 #undef __FUNC__ 807 #define __FUNC__ "MatMultTransAdd_SeqAIJ" 808 int MatMultTransAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 809 { 810 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 811 Scalar *x, *y, *v, alpha; 812 int m = a->m, n, i, *idx,shift = a->indexshift; 813 814 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 815 if (zz != yy) VecCopy(zz,yy); 816 y = y + shift; /* shift for Fortran start by 1 indexing */ 817 for ( i=0; i<m; i++ ) { 818 idx = a->j + a->i[i] + shift; 819 v = a->a + a->i[i] + shift; 820 n = a->i[i+1] - a->i[i]; 821 alpha = x[i]; 822 while (n-->0) {y[*idx++] += alpha * *v++;} 823 } 824 PLogFlops(2*a->nz); 825 return 0; 826 } 827 828 #undef __FUNC__ 829 #define __FUNC__ "MatMult_SeqAIJ" 830 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 831 { 832 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 833 Scalar *x, *y, *v, sum; 834 int m = a->m, n, i, *idx, shift = a->indexshift,*ii,jrow,j; 835 836 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 837 x = x + shift; /* shift for Fortran start by 1 indexing */ 838 idx = a->j; 839 v = a->a + shift; /* shift for Fortran start by 1 indexing */ 840 ii = a->i; 841 #if defined(USE_FORTRAN_KERNELS) 842 fortranmultaij_(&m,x,ii,idx,v,y); 843 #else 844 for ( i=0; i<m; i++ ) { 845 jrow = ii[i]; 846 n = ii[i+1] - jrow; 847 sum = 0.0; 848 for ( j=0; j<n; j++) { 849 sum += v[jrow]*x[idx[jrow]]; jrow++; 850 } 851 y[i] = sum; 852 } 853 #endif 854 PLogFlops(2*a->nz - m); 855 return 0; 856 } 857 858 #undef __FUNC__ 859 #define __FUNC__ "MatMultAdd_SeqAIJ" 860 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 861 { 862 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 863 Scalar *x, *y, *z, *v, sum; 864 int m = a->m, n, i, *idx, shift = a->indexshift,*ii,jrow,j; 865 866 867 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); VecGetArray_Fast(zz,z); 868 x = x + shift; /* shift for Fortran start by 1 indexing */ 869 idx = a->j; 870 v = a->a + shift; /* shift for Fortran start by 1 indexing */ 871 ii = a->i; 872 for ( i=0; i<m; i++ ) { 873 jrow = ii[i]; 874 n = ii[i+1] - jrow; 875 sum = y[i]; 876 for ( j=0; j<n; j++) { 877 sum += v[jrow]*x[idx[jrow]]; jrow++; 878 } 879 z[i] = sum; 880 } 881 PLogFlops(2*a->nz); 882 return 0; 883 } 884 885 /* 886 Adds diagonal pointers to sparse matrix structure. 887 */ 888 889 #undef __FUNC__ 890 #define __FUNC__ "MatMarkDiag_SeqAIJ" 891 int MatMarkDiag_SeqAIJ(Mat A) 892 { 893 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 894 int i,j, *diag, m = a->m,shift = a->indexshift; 895 896 diag = (int *) PetscMalloc( (m+1)*sizeof(int)); CHKPTRQ(diag); 897 PLogObjectMemory(A,(m+1)*sizeof(int)); 898 for ( i=0; i<a->m; i++ ) { 899 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 900 if (a->j[j]+shift == i) { 901 diag[i] = j - shift; 902 break; 903 } 904 } 905 } 906 a->diag = diag; 907 return 0; 908 } 909 910 #undef __FUNC__ 911 #define __FUNC__ "MatRelax_SeqAIJ" 912 int MatRelax_SeqAIJ(Mat A,Vec bb,double omega,MatSORType flag, 913 double fshift,int its,Vec xx) 914 { 915 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 916 Scalar *x, *b, *bs, d, *xs, sum, *v = a->a,*t,scale,*ts, *xb; 917 int ierr, *idx, *diag,n = a->n, m = a->m, i, shift = a->indexshift; 918 919 VecGetArray_Fast(xx,x); VecGetArray_Fast(bb,b); 920 if (!a->diag) {ierr = MatMarkDiag_SeqAIJ(A);CHKERRQ(ierr);} 921 diag = a->diag; 922 xs = x + shift; /* shifted by one for index start of a or a->j*/ 923 if (flag == SOR_APPLY_UPPER) { 924 /* apply ( U + D/omega) to the vector */ 925 bs = b + shift; 926 for ( i=0; i<m; i++ ) { 927 d = fshift + a->a[diag[i] + shift]; 928 n = a->i[i+1] - diag[i] - 1; 929 idx = a->j + diag[i] + (!shift); 930 v = a->a + diag[i] + (!shift); 931 sum = b[i]*d/omega; 932 SPARSEDENSEDOT(sum,bs,v,idx,n); 933 x[i] = sum; 934 } 935 return 0; 936 } 937 if (flag == SOR_APPLY_LOWER) { 938 SETERRQ(1,0,"SOR_APPLY_LOWER is not done"); 939 } 940 else if (flag & SOR_EISENSTAT) { 941 /* Let A = L + U + D; where L is lower trianglar, 942 U is upper triangular, E is diagonal; This routine applies 943 944 (L + E)^{-1} A (U + E)^{-1} 945 946 to a vector efficiently using Eisenstat's trick. This is for 947 the case of SSOR preconditioner, so E is D/omega where omega 948 is the relaxation factor. 949 */ 950 t = (Scalar *) PetscMalloc( m*sizeof(Scalar) ); CHKPTRQ(t); 951 scale = (2.0/omega) - 1.0; 952 953 /* x = (E + U)^{-1} b */ 954 for ( i=m-1; i>=0; i-- ) { 955 d = fshift + a->a[diag[i] + shift]; 956 n = a->i[i+1] - diag[i] - 1; 957 idx = a->j + diag[i] + (!shift); 958 v = a->a + diag[i] + (!shift); 959 sum = b[i]; 960 SPARSEDENSEMDOT(sum,xs,v,idx,n); 961 x[i] = omega*(sum/d); 962 } 963 964 /* t = b - (2*E - D)x */ 965 v = a->a; 966 for ( i=0; i<m; i++ ) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 967 968 /* t = (E + L)^{-1}t */ 969 ts = t + shift; /* shifted by one for index start of a or a->j*/ 970 diag = a->diag; 971 for ( i=0; i<m; i++ ) { 972 d = fshift + a->a[diag[i]+shift]; 973 n = diag[i] - a->i[i]; 974 idx = a->j + a->i[i] + shift; 975 v = a->a + a->i[i] + shift; 976 sum = t[i]; 977 SPARSEDENSEMDOT(sum,ts,v,idx,n); 978 t[i] = omega*(sum/d); 979 } 980 981 /* x = x + t */ 982 for ( i=0; i<m; i++ ) { x[i] += t[i]; } 983 PetscFree(t); 984 return 0; 985 } 986 if (flag & SOR_ZERO_INITIAL_GUESS) { 987 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 988 for ( i=0; i<m; i++ ) { 989 d = fshift + a->a[diag[i]+shift]; 990 n = diag[i] - a->i[i]; 991 idx = a->j + a->i[i] + shift; 992 v = a->a + a->i[i] + shift; 993 sum = b[i]; 994 SPARSEDENSEMDOT(sum,xs,v,idx,n); 995 x[i] = omega*(sum/d); 996 } 997 xb = x; 998 } 999 else xb = b; 1000 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1001 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1002 for ( i=0; i<m; i++ ) { 1003 x[i] *= a->a[diag[i]+shift]; 1004 } 1005 } 1006 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1007 for ( i=m-1; i>=0; i-- ) { 1008 d = fshift + a->a[diag[i] + shift]; 1009 n = a->i[i+1] - diag[i] - 1; 1010 idx = a->j + diag[i] + (!shift); 1011 v = a->a + diag[i] + (!shift); 1012 sum = xb[i]; 1013 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1014 x[i] = omega*(sum/d); 1015 } 1016 } 1017 its--; 1018 } 1019 while (its--) { 1020 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1021 for ( i=0; i<m; i++ ) { 1022 d = fshift + a->a[diag[i]+shift]; 1023 n = a->i[i+1] - a->i[i]; 1024 idx = a->j + a->i[i] + shift; 1025 v = a->a + a->i[i] + shift; 1026 sum = b[i]; 1027 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1028 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1029 } 1030 } 1031 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1032 for ( i=m-1; i>=0; i-- ) { 1033 d = fshift + a->a[diag[i] + shift]; 1034 n = a->i[i+1] - a->i[i]; 1035 idx = a->j + a->i[i] + shift; 1036 v = a->a + a->i[i] + shift; 1037 sum = b[i]; 1038 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1039 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1040 } 1041 } 1042 } 1043 return 0; 1044 } 1045 1046 #undef __FUNC__ 1047 #define __FUNC__ "MatGetInfo_SeqAIJ" 1048 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1049 { 1050 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1051 1052 info->rows_global = (double)a->m; 1053 info->columns_global = (double)a->n; 1054 info->rows_local = (double)a->m; 1055 info->columns_local = (double)a->n; 1056 info->block_size = 1.0; 1057 info->nz_allocated = (double)a->maxnz; 1058 info->nz_used = (double)a->nz; 1059 info->nz_unneeded = (double)(a->maxnz - a->nz); 1060 /* if (info->nz_unneeded != A->info.nz_unneeded) 1061 printf("space descrepancy: maxnz-nz = %d, nz_unneeded = %d\n",(int)info->nz_unneeded,(int)A->info.nz_unneeded); */ 1062 info->assemblies = (double)A->num_ass; 1063 info->mallocs = (double)a->reallocs; 1064 info->memory = A->mem; 1065 if (A->factor) { 1066 info->fill_ratio_given = A->info.fill_ratio_given; 1067 info->fill_ratio_needed = A->info.fill_ratio_needed; 1068 info->factor_mallocs = A->info.factor_mallocs; 1069 } else { 1070 info->fill_ratio_given = 0; 1071 info->fill_ratio_needed = 0; 1072 info->factor_mallocs = 0; 1073 } 1074 return 0; 1075 } 1076 1077 extern int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,double,Mat*); 1078 extern int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1079 extern int MatLUFactor_SeqAIJ(Mat,IS,IS,double); 1080 extern int MatSolve_SeqAIJ(Mat,Vec,Vec); 1081 extern int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1082 extern int MatSolveTrans_SeqAIJ(Mat,Vec,Vec); 1083 extern int MatSolveTransAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1084 1085 #undef __FUNC__ 1086 #define __FUNC__ "MatZeroRows_SeqAIJ" 1087 int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 1088 { 1089 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1090 int i,ierr,N, *rows,m = a->m - 1,shift = a->indexshift; 1091 1092 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 1093 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 1094 if (diag) { 1095 for ( i=0; i<N; i++ ) { 1096 if (rows[i] < 0 || rows[i] > m) SETERRQ(1,0,"row out of range"); 1097 if (a->ilen[rows[i]] > 0) { /* in case row was completely empty */ 1098 a->ilen[rows[i]] = 1; 1099 a->a[a->i[rows[i]]+shift] = *diag; 1100 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1101 } 1102 else { 1103 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES); 1104 CHKERRQ(ierr); 1105 } 1106 } 1107 } 1108 else { 1109 for ( i=0; i<N; i++ ) { 1110 if (rows[i] < 0 || rows[i] > m) SETERRQ(1,0,"row out of range"); 1111 a->ilen[rows[i]] = 0; 1112 } 1113 } 1114 ISRestoreIndices(is,&rows); 1115 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1116 return 0; 1117 } 1118 1119 #undef __FUNC__ 1120 #define __FUNC__ "MatGetSize_SeqAIJ" 1121 int MatGetSize_SeqAIJ(Mat A,int *m,int *n) 1122 { 1123 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1124 *m = a->m; *n = a->n; 1125 return 0; 1126 } 1127 1128 #undef __FUNC__ 1129 #define __FUNC__ "MatGetOwnershipRange_SeqAIJ" 1130 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1131 { 1132 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1133 *m = 0; *n = a->m; 1134 return 0; 1135 } 1136 1137 #undef __FUNC__ 1138 #define __FUNC__ "MatGetRow_SeqAIJ" 1139 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1140 { 1141 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1142 int *itmp,i,shift = a->indexshift; 1143 1144 if (row < 0 || row >= a->m) SETERRQ(1,0,"Row out of range"); 1145 1146 *nz = a->i[row+1] - a->i[row]; 1147 if (v) *v = a->a + a->i[row] + shift; 1148 if (idx) { 1149 itmp = a->j + a->i[row] + shift; 1150 if (*nz && shift) { 1151 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 1152 for ( i=0; i<(*nz); i++ ) {(*idx)[i] = itmp[i] + shift;} 1153 } else if (*nz) { 1154 *idx = itmp; 1155 } 1156 else *idx = 0; 1157 } 1158 return 0; 1159 } 1160 1161 #undef __FUNC__ 1162 #define __FUNC__ "MatRestoreRow_SeqAIJ" 1163 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1164 { 1165 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1166 if (idx) {if (*idx && a->indexshift) PetscFree(*idx);} 1167 return 0; 1168 } 1169 1170 #undef __FUNC__ 1171 #define __FUNC__ "MatNorm_SeqAIJ" 1172 int MatNorm_SeqAIJ(Mat A,NormType type,double *norm) 1173 { 1174 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1175 Scalar *v = a->a; 1176 double sum = 0.0; 1177 int i, j,shift = a->indexshift; 1178 1179 if (type == NORM_FROBENIUS) { 1180 for (i=0; i<a->nz; i++ ) { 1181 #if defined(PETSC_COMPLEX) 1182 sum += real(conj(*v)*(*v)); v++; 1183 #else 1184 sum += (*v)*(*v); v++; 1185 #endif 1186 } 1187 *norm = sqrt(sum); 1188 } 1189 else if (type == NORM_1) { 1190 double *tmp; 1191 int *jj = a->j; 1192 tmp = (double *) PetscMalloc( (a->n+1)*sizeof(double) ); CHKPTRQ(tmp); 1193 PetscMemzero(tmp,a->n*sizeof(double)); 1194 *norm = 0.0; 1195 for ( j=0; j<a->nz; j++ ) { 1196 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1197 } 1198 for ( j=0; j<a->n; j++ ) { 1199 if (tmp[j] > *norm) *norm = tmp[j]; 1200 } 1201 PetscFree(tmp); 1202 } 1203 else if (type == NORM_INFINITY) { 1204 *norm = 0.0; 1205 for ( j=0; j<a->m; j++ ) { 1206 v = a->a + a->i[j] + shift; 1207 sum = 0.0; 1208 for ( i=0; i<a->i[j+1]-a->i[j]; i++ ) { 1209 sum += PetscAbsScalar(*v); v++; 1210 } 1211 if (sum > *norm) *norm = sum; 1212 } 1213 } 1214 else { 1215 SETERRQ(1,0,"No support for two norm yet"); 1216 } 1217 return 0; 1218 } 1219 1220 #undef __FUNC__ 1221 #define __FUNC__ "MatTranspose_SeqAIJ" 1222 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1223 { 1224 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1225 Mat C; 1226 int i, ierr, *aj = a->j, *ai = a->i, m = a->m, len, *col; 1227 int shift = a->indexshift; 1228 Scalar *array = a->a; 1229 1230 if (B == PETSC_NULL && m != a->n) 1231 SETERRQ(1,0,"Square matrix only for in-place"); 1232 col = (int *) PetscMalloc((1+a->n)*sizeof(int)); CHKPTRQ(col); 1233 PetscMemzero(col,(1+a->n)*sizeof(int)); 1234 if (shift) { 1235 for ( i=0; i<ai[m]-1; i++ ) aj[i] -= 1; 1236 } 1237 for ( i=0; i<ai[m]+shift; i++ ) col[aj[i]] += 1; 1238 ierr = MatCreateSeqAIJ(A->comm,a->n,m,0,col,&C); CHKERRQ(ierr); 1239 PetscFree(col); 1240 for ( i=0; i<m; i++ ) { 1241 len = ai[i+1]-ai[i]; 1242 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES); CHKERRQ(ierr); 1243 array += len; aj += len; 1244 } 1245 if (shift) { 1246 for ( i=0; i<ai[m]-1; i++ ) aj[i] += 1; 1247 } 1248 1249 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1250 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1251 1252 if (B != PETSC_NULL) { 1253 *B = C; 1254 } else { 1255 /* This isn't really an in-place transpose */ 1256 PetscFree(a->a); 1257 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 1258 if (a->diag) PetscFree(a->diag); 1259 if (a->ilen) PetscFree(a->ilen); 1260 if (a->imax) PetscFree(a->imax); 1261 if (a->solve_work) PetscFree(a->solve_work); 1262 if (a->inode.size) PetscFree(a->inode.size); 1263 PetscFree(a); 1264 PetscMemcpy(A,C,sizeof(struct _p_Mat)); 1265 PetscHeaderDestroy(C); 1266 } 1267 return 0; 1268 } 1269 1270 #undef __FUNC__ 1271 #define __FUNC__ "MatDiagonalScale_SeqAIJ" 1272 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1273 { 1274 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1275 Scalar *l,*r,x,*v; 1276 int i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1277 1278 if (ll) { 1279 /* The local size is used so that VecMPI can be passed to this routine 1280 by MatDiagonalScale_MPIAIJ */ 1281 VecGetLocalSize_Fast(ll,m); 1282 if (m != a->m) SETERRQ(1,0,"Left scaling vector wrong length"); 1283 VecGetArray_Fast(ll,l); 1284 v = a->a; 1285 for ( i=0; i<m; i++ ) { 1286 x = l[i]; 1287 M = a->i[i+1] - a->i[i]; 1288 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1289 } 1290 PLogFlops(nz); 1291 } 1292 if (rr) { 1293 VecGetLocalSize_Fast(rr,n); 1294 if (n != a->n) SETERRQ(1,0,"Right scaling vector wrong length"); 1295 VecGetArray_Fast(rr,r); 1296 v = a->a; jj = a->j; 1297 for ( i=0; i<nz; i++ ) { 1298 (*v++) *= r[*jj++ + shift]; 1299 } 1300 PLogFlops(nz); 1301 } 1302 return 0; 1303 } 1304 1305 #undef __FUNC__ 1306 #define __FUNC__ "MatGetSubMatrix_SeqAIJ" 1307 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,MatGetSubMatrixCall scall,Mat *B) 1308 { 1309 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1310 int nznew, *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1311 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1312 register int sum,lensi; 1313 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1314 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1315 Scalar *a_new,*mat_a; 1316 Mat C; 1317 1318 ierr = ISSorted(isrow,(PetscTruth*)&i); 1319 if (!i) SETERRQ(1,0,"ISrow is not sorted"); 1320 ierr = ISSorted(iscol,(PetscTruth*)&i); 1321 if (!i) SETERRQ(1,0,"IScol is not sorted"); 1322 1323 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1324 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1325 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1326 1327 if (ISStrideGetInfo(iscol,&first,&step) && step == 1) { /* no need to sort */ 1328 /* special case of contiguous rows */ 1329 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1330 starts = lens + ncols; 1331 /* loop over new rows determining lens and starting points */ 1332 for (i=0; i<nrows; i++) { 1333 kstart = ai[irow[i]]+shift; 1334 kend = kstart + ailen[irow[i]]; 1335 for ( k=kstart; k<kend; k++ ) { 1336 if (aj[k]+shift >= first) { 1337 starts[i] = k; 1338 break; 1339 } 1340 } 1341 sum = 0; 1342 while (k < kend) { 1343 if (aj[k++]+shift >= first+ncols) break; 1344 sum++; 1345 } 1346 lens[i] = sum; 1347 } 1348 /* create submatrix */ 1349 if (scall == MAT_REUSE_MATRIX) { 1350 int n_cols,n_rows; 1351 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1352 if (n_rows != nrows || n_cols != ncols) SETERRQ(1,0,""); 1353 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1354 C = *B; 1355 } 1356 else { 1357 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1358 } 1359 c = (Mat_SeqAIJ*) C->data; 1360 1361 /* loop over rows inserting into submatrix */ 1362 a_new = c->a; 1363 j_new = c->j; 1364 i_new = c->i; 1365 i_new[0] = -shift; 1366 for (i=0; i<nrows; i++) { 1367 ii = starts[i]; 1368 lensi = lens[i]; 1369 for ( k=0; k<lensi; k++ ) { 1370 *j_new++ = aj[ii+k] - first; 1371 } 1372 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1373 a_new += lensi; 1374 i_new[i+1] = i_new[i] + lensi; 1375 c->ilen[i] = lensi; 1376 } 1377 PetscFree(lens); 1378 } 1379 else { 1380 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1381 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1382 ssmap = smap + shift; 1383 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1384 PetscMemzero(smap,oldcols*sizeof(int)); 1385 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1386 /* determine lens of each row */ 1387 for (i=0; i<nrows; i++) { 1388 kstart = ai[irow[i]]+shift; 1389 kend = kstart + a->ilen[irow[i]]; 1390 lens[i] = 0; 1391 for ( k=kstart; k<kend; k++ ) { 1392 if (ssmap[aj[k]]) { 1393 lens[i]++; 1394 } 1395 } 1396 } 1397 /* Create and fill new matrix */ 1398 if (scall == MAT_REUSE_MATRIX) { 1399 c = (Mat_SeqAIJ *)((*B)->data); 1400 1401 if (c->m != nrows || c->n != ncols) SETERRQ(1,0,""); 1402 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1403 SETERRQ(1,0,"Cannot reuse matrix. wrong no of nonzeros"); 1404 } 1405 PetscMemzero(c->ilen,c->m*sizeof(int)); 1406 C = *B; 1407 } 1408 else { 1409 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1410 } 1411 c = (Mat_SeqAIJ *)(C->data); 1412 for (i=0; i<nrows; i++) { 1413 row = irow[i]; 1414 nznew = 0; 1415 kstart = ai[row]+shift; 1416 kend = kstart + a->ilen[row]; 1417 mat_i = c->i[i]+shift; 1418 mat_j = c->j + mat_i; 1419 mat_a = c->a + mat_i; 1420 mat_ilen = c->ilen + i; 1421 for ( k=kstart; k<kend; k++ ) { 1422 if ((tcol=ssmap[a->j[k]])) { 1423 *mat_j++ = tcol - (!shift); 1424 *mat_a++ = a->a[k]; 1425 (*mat_ilen)++; 1426 1427 } 1428 } 1429 } 1430 /* Free work space */ 1431 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1432 PetscFree(smap); PetscFree(lens); 1433 } 1434 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1435 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1436 1437 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1438 *B = C; 1439 return 0; 1440 } 1441 1442 /* 1443 note: This can only work for identity for row and col. It would 1444 be good to check this and otherwise generate an error. 1445 */ 1446 #undef __FUNC__ 1447 #define __FUNC__ "MatILUFactor_SeqAIJ" 1448 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1449 { 1450 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1451 int ierr; 1452 Mat outA; 1453 1454 if (fill != 0) SETERRQ(1,0,"Only fill=0 supported"); 1455 1456 outA = inA; 1457 inA->factor = FACTOR_LU; 1458 a->row = row; 1459 a->col = col; 1460 1461 if (!a->solve_work) { /* this matrix may have been factored before */ 1462 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar)); CHKPTRQ(a->solve_work); 1463 } 1464 1465 if (!a->diag) { 1466 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1467 } 1468 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1469 return 0; 1470 } 1471 1472 #include "pinclude/plapack.h" 1473 #undef __FUNC__ 1474 #define __FUNC__ "MatScale_SeqAIJ" 1475 int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1476 { 1477 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1478 int one = 1; 1479 BLscal_( &a->nz, alpha, a->a, &one ); 1480 PLogFlops(a->nz); 1481 return 0; 1482 } 1483 1484 #undef __FUNC__ 1485 #define __FUNC__ "MatGetSubMatrices_SeqAIJ" 1486 int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall, 1487 Mat **B) 1488 { 1489 int ierr,i; 1490 1491 if (scall == MAT_INITIAL_MATRIX) { 1492 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1493 } 1494 1495 for ( i=0; i<n; i++ ) { 1496 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 1497 } 1498 return 0; 1499 } 1500 1501 #undef __FUNC__ 1502 #define __FUNC__ "MatGetBlockSize_SeqAIJ" 1503 int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1504 { 1505 *bs = 1; 1506 return 0; 1507 } 1508 1509 #undef __FUNC__ 1510 #define __FUNC__ "MatIncreaseOverlap_SeqAIJ" 1511 int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1512 { 1513 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1514 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1515 int start, end, *ai, *aj; 1516 char *table; 1517 shift = a->indexshift; 1518 m = a->m; 1519 ai = a->i; 1520 aj = a->j+shift; 1521 1522 if (ov < 0) SETERRQ(1,0,"illegal overlap value used"); 1523 1524 table = (char *) PetscMalloc((m/BITSPERBYTE +1)*sizeof(char)); CHKPTRQ(table); 1525 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1526 1527 for ( i=0; i<is_max; i++ ) { 1528 /* Initialize the two local arrays */ 1529 isz = 0; 1530 PetscMemzero(table,(m/BITSPERBYTE +1)*sizeof(char)); 1531 1532 /* Extract the indices, assume there can be duplicate entries */ 1533 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1534 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1535 1536 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1537 for ( j=0; j<n ; ++j){ 1538 if(!BT_LOOKUP(table, idx[j])) { nidx[isz++] = idx[j];} 1539 } 1540 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1541 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1542 1543 k = 0; 1544 for ( j=0; j<ov; j++){ /* for each overlap */ 1545 n = isz; 1546 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1547 row = nidx[k]; 1548 start = ai[row]; 1549 end = ai[row+1]; 1550 for ( l = start; l<end ; l++){ 1551 val = aj[l] + shift; 1552 if (!BT_LOOKUP(table,val)) {nidx[isz++] = val;} 1553 } 1554 } 1555 } 1556 ierr = ISCreateGeneral(PETSC_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1557 } 1558 PetscFree(table); 1559 PetscFree(nidx); 1560 return 0; 1561 } 1562 1563 /* -------------------------------------------------------------- */ 1564 #undef __FUNC__ 1565 #define __FUNC__ "MatPermute_SeqAIJ" 1566 int MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 1567 { 1568 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1569 Scalar *vwork; 1570 int i, ierr, nz, m = a->m, n = a->n, *cwork; 1571 int *row,*col,*cnew,j,*lens; 1572 IS icolp,irowp; 1573 1574 ierr = ISInvertPermutation(rowp,&irowp); CHKERRQ(ierr); 1575 ierr = ISGetIndices(irowp,&row); CHKERRQ(ierr); 1576 ierr = ISInvertPermutation(colp,&icolp); CHKERRQ(ierr); 1577 ierr = ISGetIndices(icolp,&col); CHKERRQ(ierr); 1578 1579 /* determine lengths of permuted rows */ 1580 lens = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(lens); 1581 for (i=0; i<m; i++ ) { 1582 lens[row[i]] = a->i[i+1] - a->i[i]; 1583 } 1584 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1585 PetscFree(lens); 1586 1587 cnew = (int *) PetscMalloc( n*sizeof(int) ); CHKPTRQ(cnew); 1588 for (i=0; i<m; i++) { 1589 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1590 for (j=0; j<nz; j++ ) { cnew[j] = col[cwork[j]];} 1591 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES); CHKERRQ(ierr); 1592 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1593 } 1594 PetscFree(cnew); 1595 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1596 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1597 ierr = ISRestoreIndices(irowp,&row); CHKERRQ(ierr); 1598 ierr = ISRestoreIndices(icolp,&col); CHKERRQ(ierr); 1599 ierr = ISDestroy(irowp); CHKERRQ(ierr); 1600 ierr = ISDestroy(icolp); CHKERRQ(ierr); 1601 return 0; 1602 } 1603 1604 #undef __FUNC__ 1605 #define __FUNC__ "MatPrintHelp_SeqAIJ" 1606 int MatPrintHelp_SeqAIJ(Mat A) 1607 { 1608 static int called = 0; 1609 MPI_Comm comm = A->comm; 1610 1611 if (called) return 0; else called = 1; 1612 PetscPrintf(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1613 PetscPrintf(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n"); 1614 PetscPrintf(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n"); 1615 PetscPrintf(comm," -mat_aij_no_inode: Do not use inodes\n"); 1616 PetscPrintf(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n"); 1617 #if defined(HAVE_ESSL) 1618 PetscPrintf(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n"); 1619 #endif 1620 return 0; 1621 } 1622 extern int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1623 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1624 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1625 1626 /* -------------------------------------------------------------------*/ 1627 static struct _MatOps MatOps = {MatSetValues_SeqAIJ, 1628 MatGetRow_SeqAIJ,MatRestoreRow_SeqAIJ, 1629 MatMult_SeqAIJ,MatMultAdd_SeqAIJ, 1630 MatMultTrans_SeqAIJ,MatMultTransAdd_SeqAIJ, 1631 MatSolve_SeqAIJ,MatSolveAdd_SeqAIJ, 1632 MatSolveTrans_SeqAIJ,MatSolveTransAdd_SeqAIJ, 1633 MatLUFactor_SeqAIJ,0, 1634 MatRelax_SeqAIJ, 1635 MatTranspose_SeqAIJ, 1636 MatGetInfo_SeqAIJ,MatEqual_SeqAIJ, 1637 MatGetDiagonal_SeqAIJ,MatDiagonalScale_SeqAIJ,MatNorm_SeqAIJ, 1638 0,MatAssemblyEnd_SeqAIJ, 1639 MatCompress_SeqAIJ, 1640 MatSetOption_SeqAIJ,MatZeroEntries_SeqAIJ,MatZeroRows_SeqAIJ, 1641 MatLUFactorSymbolic_SeqAIJ,MatLUFactorNumeric_SeqAIJ,0,0, 1642 MatGetSize_SeqAIJ,MatGetSize_SeqAIJ,MatGetOwnershipRange_SeqAIJ, 1643 MatILUFactorSymbolic_SeqAIJ,0, 1644 0,0, 1645 MatConvertSameType_SeqAIJ,0,0, 1646 MatILUFactor_SeqAIJ,0,0, 1647 MatGetSubMatrices_SeqAIJ,MatIncreaseOverlap_SeqAIJ, 1648 MatGetValues_SeqAIJ,0, 1649 MatPrintHelp_SeqAIJ, 1650 MatScale_SeqAIJ,0,0, 1651 MatILUDTFactor_SeqAIJ, 1652 MatGetBlockSize_SeqAIJ, 1653 MatGetRowIJ_SeqAIJ, 1654 MatRestoreRowIJ_SeqAIJ, 1655 MatGetColumnIJ_SeqAIJ, 1656 MatRestoreColumnIJ_SeqAIJ, 1657 MatFDColoringCreate_SeqAIJ, 1658 MatColoringPatch_SeqAIJ, 1659 0, 1660 MatPermute_SeqAIJ}; 1661 1662 extern int MatUseSuperLU_SeqAIJ(Mat); 1663 extern int MatUseEssl_SeqAIJ(Mat); 1664 extern int MatUseDXML_SeqAIJ(Mat); 1665 1666 #undef __FUNC__ 1667 #define __FUNC__ "MatCreateSeqAIJ" 1668 /*@C 1669 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 1670 (the default parallel PETSc format). For good matrix assembly performance 1671 the user should preallocate the matrix storage by setting the parameter nz 1672 (or the array nzz). By setting these parameters accurately, performance 1673 during matrix assembly can be increased by more than a factor of 50. 1674 1675 Input Parameters: 1676 . comm - MPI communicator, set to PETSC_COMM_SELF 1677 . m - number of rows 1678 . n - number of columns 1679 . nz - number of nonzeros per row (same for all rows) 1680 . nzz - array containing the number of nonzeros in the various rows 1681 (possibly different for each row) or PETSC_NULL 1682 1683 Output Parameter: 1684 . A - the matrix 1685 1686 Notes: 1687 The AIJ format (also called the Yale sparse matrix format or 1688 compressed row storage), is fully compatible with standard Fortran 77 1689 storage. That is, the stored row and column indices can begin at 1690 either one (as in Fortran) or zero. See the users' manual for details. 1691 1692 Specify the preallocated storage with either nz or nnz (not both). 1693 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1694 allocation. For large problems you MUST preallocate memory or you 1695 will get TERRIBLE performance, see the users' manual chapter on matrices. 1696 1697 By default, this format uses inodes (identical nodes) when possible, to 1698 improve numerical efficiency of matrix-vector products and solves. We 1699 search for consecutive rows with the same nonzero structure, thereby 1700 reusing matrix information to achieve increased efficiency. 1701 1702 Options Database Keys: 1703 $ -mat_aij_no_inode - Do not use inodes 1704 $ -mat_aij_inode_limit <limit> - Set inode limit. 1705 $ (max limit=5) 1706 $ -mat_aij_oneindex - Internally use indexing starting at 1 1707 $ rather than 0. Note: When calling MatSetValues(), 1708 $ the user still MUST index entries starting at 0! 1709 1710 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues() 1711 @*/ 1712 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 1713 { 1714 Mat B; 1715 Mat_SeqAIJ *b; 1716 int i, len, ierr, flg,size; 1717 1718 MPI_Comm_size(comm,&size); 1719 if (size > 1) SETERRQ(1,0,"Comm must be of size 1"); 1720 1721 *A = 0; 1722 PetscHeaderCreate(B,_p_Mat,MAT_COOKIE,MATSEQAIJ,comm,MatDestroy,MatView); 1723 PLogObjectCreate(B); 1724 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 1725 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 1726 PetscMemcpy(&B->ops,&MatOps,sizeof(struct _MatOps)); 1727 B->destroy = MatDestroy_SeqAIJ; 1728 B->view = MatView_SeqAIJ; 1729 B->factor = 0; 1730 B->lupivotthreshold = 1.0; 1731 B->mapping = 0; 1732 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold, 1733 &flg); CHKERRQ(ierr); 1734 b->ilu_preserve_row_sums = PETSC_FALSE; 1735 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums", 1736 (int*) &b->ilu_preserve_row_sums); CHKERRQ(ierr); 1737 b->row = 0; 1738 b->col = 0; 1739 b->indexshift = 0; 1740 b->reallocs = 0; 1741 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 1742 if (flg) b->indexshift = -1; 1743 1744 b->m = m; B->m = m; B->M = m; 1745 b->n = n; B->n = n; B->N = n; 1746 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 1747 if (nnz == PETSC_NULL) { 1748 if (nz == PETSC_DEFAULT) nz = 10; 1749 else if (nz <= 0) nz = 1; 1750 for ( i=0; i<m; i++ ) b->imax[i] = nz; 1751 nz = nz*m; 1752 } 1753 else { 1754 nz = 0; 1755 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1756 } 1757 1758 /* allocate the matrix space */ 1759 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 1760 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1761 b->j = (int *) (b->a + nz); 1762 PetscMemzero(b->j,nz*sizeof(int)); 1763 b->i = b->j + nz; 1764 b->singlemalloc = 1; 1765 1766 b->i[0] = -b->indexshift; 1767 for (i=1; i<m+1; i++) { 1768 b->i[i] = b->i[i-1] + b->imax[i-1]; 1769 } 1770 1771 /* b->ilen will count nonzeros in each row so far. */ 1772 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 1773 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1774 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 1775 1776 b->nz = 0; 1777 b->maxnz = nz; 1778 b->sorted = 0; 1779 b->roworiented = 1; 1780 b->nonew = 0; 1781 b->diag = 0; 1782 b->solve_work = 0; 1783 b->spptr = 0; 1784 b->inode.node_count = 0; 1785 b->inode.size = 0; 1786 b->inode.limit = 5; 1787 b->inode.max_limit = 5; 1788 B->info.nz_unneeded = (double)b->maxnz; 1789 1790 *A = B; 1791 1792 /* SuperLU is not currently supported through PETSc */ 1793 #if defined(HAVE_SUPERLU) 1794 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 1795 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 1796 #endif 1797 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 1798 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 1799 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 1800 if (flg) { 1801 if (!b->indexshift) SETERRQ(1,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 1802 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 1803 } 1804 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1805 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1806 return 0; 1807 } 1808 1809 #undef __FUNC__ 1810 #define __FUNC__ "MatConvertSameType_SeqAIJ" 1811 int MatConvertSameType_SeqAIJ(Mat A,Mat *B,int cpvalues) 1812 { 1813 Mat C; 1814 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 1815 int i,len, m = a->m,shift = a->indexshift; 1816 1817 *B = 0; 1818 PetscHeaderCreate(C,_p_Mat,MAT_COOKIE,MATSEQAIJ,A->comm,MatDestroy,MatView); 1819 PLogObjectCreate(C); 1820 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 1821 PetscMemcpy(&C->ops,&A->ops,sizeof(struct _MatOps)); 1822 C->destroy = MatDestroy_SeqAIJ; 1823 C->view = MatView_SeqAIJ; 1824 C->factor = A->factor; 1825 c->row = 0; 1826 c->col = 0; 1827 c->indexshift = shift; 1828 C->assembled = PETSC_TRUE; 1829 1830 c->m = C->m = a->m; 1831 c->n = C->n = a->n; 1832 C->M = a->m; 1833 C->N = a->n; 1834 1835 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 1836 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 1837 for ( i=0; i<m; i++ ) { 1838 c->imax[i] = a->imax[i]; 1839 c->ilen[i] = a->ilen[i]; 1840 } 1841 1842 /* allocate the matrix space */ 1843 c->singlemalloc = 1; 1844 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 1845 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1846 c->j = (int *) (c->a + a->i[m] + shift); 1847 c->i = c->j + a->i[m] + shift; 1848 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 1849 if (m > 0) { 1850 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 1851 if (cpvalues == COPY_VALUES) { 1852 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 1853 } 1854 } 1855 1856 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1857 c->sorted = a->sorted; 1858 c->roworiented = a->roworiented; 1859 c->nonew = a->nonew; 1860 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 1861 1862 if (a->diag) { 1863 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 1864 PLogObjectMemory(C,(m+1)*sizeof(int)); 1865 for ( i=0; i<m; i++ ) { 1866 c->diag[i] = a->diag[i]; 1867 } 1868 } 1869 else c->diag = 0; 1870 c->inode.limit = a->inode.limit; 1871 c->inode.max_limit = a->inode.max_limit; 1872 if (a->inode.size){ 1873 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 1874 c->inode.node_count = a->inode.node_count; 1875 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 1876 } else { 1877 c->inode.size = 0; 1878 c->inode.node_count = 0; 1879 } 1880 c->nz = a->nz; 1881 c->maxnz = a->maxnz; 1882 c->solve_work = 0; 1883 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1884 1885 *B = C; 1886 return 0; 1887 } 1888 1889 #undef __FUNC__ 1890 #define __FUNC__ "MatLoad_SeqAIJ" 1891 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 1892 { 1893 Mat_SeqAIJ *a; 1894 Mat B; 1895 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 1896 MPI_Comm comm; 1897 1898 PetscObjectGetComm((PetscObject) viewer,&comm); 1899 MPI_Comm_size(comm,&size); 1900 if (size > 1) SETERRQ(1,0,"view must have one processor"); 1901 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 1902 ierr = PetscBinaryRead(fd,header,4,BINARY_INT); CHKERRQ(ierr); 1903 if (header[0] != MAT_COOKIE) SETERRQ(1,0,"not matrix object in file"); 1904 M = header[1]; N = header[2]; nz = header[3]; 1905 1906 /* read in row lengths */ 1907 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 1908 ierr = PetscBinaryRead(fd,rowlengths,M,BINARY_INT); CHKERRQ(ierr); 1909 1910 /* create our matrix */ 1911 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 1912 B = *A; 1913 a = (Mat_SeqAIJ *) B->data; 1914 shift = a->indexshift; 1915 1916 /* read in column indices and adjust for Fortran indexing*/ 1917 ierr = PetscBinaryRead(fd,a->j,nz,BINARY_INT); CHKERRQ(ierr); 1918 if (shift) { 1919 for ( i=0; i<nz; i++ ) { 1920 a->j[i] += 1; 1921 } 1922 } 1923 1924 /* read in nonzero values */ 1925 ierr = PetscBinaryRead(fd,a->a,nz,BINARY_SCALAR); CHKERRQ(ierr); 1926 1927 /* set matrix "i" values */ 1928 a->i[0] = -shift; 1929 for ( i=1; i<= M; i++ ) { 1930 a->i[i] = a->i[i-1] + rowlengths[i-1]; 1931 a->ilen[i-1] = rowlengths[i-1]; 1932 } 1933 PetscFree(rowlengths); 1934 1935 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1936 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1937 return 0; 1938 } 1939 1940 #undef __FUNC__ 1941 #define __FUNC__ "MatEqual_SeqAIJ" 1942 int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 1943 { 1944 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 1945 1946 if (B->type !=MATSEQAIJ)SETERRQ(1,0,"Matrices must be same type"); 1947 1948 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 1949 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 1950 (a->indexshift != b->indexshift)) { 1951 *flg = PETSC_FALSE; return 0; 1952 } 1953 1954 /* if the a->i are the same */ 1955 if (PetscMemcmp(a->i,b->i,(a->m+1)*sizeof(int))) { 1956 *flg = PETSC_FALSE; return 0; 1957 } 1958 1959 /* if a->j are the same */ 1960 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 1961 *flg = PETSC_FALSE; return 0; 1962 } 1963 1964 /* if a->a are the same */ 1965 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 1966 *flg = PETSC_FALSE; return 0; 1967 } 1968 *flg = PETSC_TRUE; 1969 return 0; 1970 1971 } 1972