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