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