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