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