1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: aij.c,v 1.236 1997/09/18 03:06:45 balay Exp bsmith $"; 3 #endif 4 5 /* 6 Defines the basic matrix operations for the AIJ (compressed row) 7 matrix storage format. 8 */ 9 10 #include "pinclude/pviewer.h" 11 #include "sys.h" 12 #include "src/mat/impls/aij/seq/aij.h" 13 #include "src/vec/vecimpl.h" 14 #include "src/inline/spops.h" 15 #include "src/inline/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,PETSC_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,PETSC_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,PETSC_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 if (m) *m = a->m; 1128 if (n) *n = a->n; 1129 return 0; 1130 } 1131 1132 #undef __FUNC__ 1133 #define __FUNC__ "MatGetOwnershipRange_SeqAIJ" 1134 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1135 { 1136 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1137 *m = 0; *n = a->m; 1138 return 0; 1139 } 1140 1141 #undef __FUNC__ 1142 #define __FUNC__ "MatGetRow_SeqAIJ" 1143 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1144 { 1145 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1146 int *itmp,i,shift = a->indexshift; 1147 1148 if (row < 0 || row >= a->m) SETERRQ(1,0,"Row out of range"); 1149 1150 *nz = a->i[row+1] - a->i[row]; 1151 if (v) *v = a->a + a->i[row] + shift; 1152 if (idx) { 1153 itmp = a->j + a->i[row] + shift; 1154 if (*nz && shift) { 1155 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 1156 for ( i=0; i<(*nz); i++ ) {(*idx)[i] = itmp[i] + shift;} 1157 } else if (*nz) { 1158 *idx = itmp; 1159 } 1160 else *idx = 0; 1161 } 1162 return 0; 1163 } 1164 1165 #undef __FUNC__ 1166 #define __FUNC__ "MatRestoreRow_SeqAIJ" 1167 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1168 { 1169 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1170 if (idx) {if (*idx && a->indexshift) PetscFree(*idx);} 1171 return 0; 1172 } 1173 1174 #undef __FUNC__ 1175 #define __FUNC__ "MatNorm_SeqAIJ" 1176 int MatNorm_SeqAIJ(Mat A,NormType type,double *norm) 1177 { 1178 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1179 Scalar *v = a->a; 1180 double sum = 0.0; 1181 int i, j,shift = a->indexshift; 1182 1183 if (type == NORM_FROBENIUS) { 1184 for (i=0; i<a->nz; i++ ) { 1185 #if defined(PETSC_COMPLEX) 1186 sum += real(conj(*v)*(*v)); v++; 1187 #else 1188 sum += (*v)*(*v); v++; 1189 #endif 1190 } 1191 *norm = sqrt(sum); 1192 } 1193 else if (type == NORM_1) { 1194 double *tmp; 1195 int *jj = a->j; 1196 tmp = (double *) PetscMalloc( (a->n+1)*sizeof(double) ); CHKPTRQ(tmp); 1197 PetscMemzero(tmp,a->n*sizeof(double)); 1198 *norm = 0.0; 1199 for ( j=0; j<a->nz; j++ ) { 1200 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1201 } 1202 for ( j=0; j<a->n; j++ ) { 1203 if (tmp[j] > *norm) *norm = tmp[j]; 1204 } 1205 PetscFree(tmp); 1206 } 1207 else if (type == NORM_INFINITY) { 1208 *norm = 0.0; 1209 for ( j=0; j<a->m; j++ ) { 1210 v = a->a + a->i[j] + shift; 1211 sum = 0.0; 1212 for ( i=0; i<a->i[j+1]-a->i[j]; i++ ) { 1213 sum += PetscAbsScalar(*v); v++; 1214 } 1215 if (sum > *norm) *norm = sum; 1216 } 1217 } 1218 else { 1219 SETERRQ(1,0,"No support for two norm yet"); 1220 } 1221 return 0; 1222 } 1223 1224 #undef __FUNC__ 1225 #define __FUNC__ "MatTranspose_SeqAIJ" 1226 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1227 { 1228 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1229 Mat C; 1230 int i, ierr, *aj = a->j, *ai = a->i, m = a->m, len, *col; 1231 int shift = a->indexshift; 1232 Scalar *array = a->a; 1233 1234 if (B == PETSC_NULL && m != a->n) 1235 SETERRQ(1,0,"Square matrix only for in-place"); 1236 col = (int *) PetscMalloc((1+a->n)*sizeof(int)); CHKPTRQ(col); 1237 PetscMemzero(col,(1+a->n)*sizeof(int)); 1238 if (shift) { 1239 for ( i=0; i<ai[m]-1; i++ ) aj[i] -= 1; 1240 } 1241 for ( i=0; i<ai[m]+shift; i++ ) col[aj[i]] += 1; 1242 ierr = MatCreateSeqAIJ(A->comm,a->n,m,0,col,&C); CHKERRQ(ierr); 1243 PetscFree(col); 1244 for ( i=0; i<m; i++ ) { 1245 len = ai[i+1]-ai[i]; 1246 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES); CHKERRQ(ierr); 1247 array += len; aj += len; 1248 } 1249 if (shift) { 1250 for ( i=0; i<ai[m]-1; i++ ) aj[i] += 1; 1251 } 1252 1253 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1254 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1255 1256 if (B != PETSC_NULL) { 1257 *B = C; 1258 } else { 1259 /* This isn't really an in-place transpose */ 1260 PetscFree(a->a); 1261 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 1262 if (a->diag) PetscFree(a->diag); 1263 if (a->ilen) PetscFree(a->ilen); 1264 if (a->imax) PetscFree(a->imax); 1265 if (a->solve_work) PetscFree(a->solve_work); 1266 if (a->inode.size) PetscFree(a->inode.size); 1267 PetscFree(a); 1268 PetscMemcpy(A,C,sizeof(struct _p_Mat)); 1269 PetscHeaderDestroy(C); 1270 } 1271 return 0; 1272 } 1273 1274 #undef __FUNC__ 1275 #define __FUNC__ "MatDiagonalScale_SeqAIJ" 1276 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1277 { 1278 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1279 Scalar *l,*r,x,*v; 1280 int i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1281 1282 if (ll) { 1283 /* The local size is used so that VecMPI can be passed to this routine 1284 by MatDiagonalScale_MPIAIJ */ 1285 VecGetLocalSize_Fast(ll,m); 1286 if (m != a->m) SETERRQ(1,0,"Left scaling vector wrong length"); 1287 VecGetArray_Fast(ll,l); 1288 v = a->a; 1289 for ( i=0; i<m; i++ ) { 1290 x = l[i]; 1291 M = a->i[i+1] - a->i[i]; 1292 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1293 } 1294 PLogFlops(nz); 1295 } 1296 if (rr) { 1297 VecGetLocalSize_Fast(rr,n); 1298 if (n != a->n) SETERRQ(1,0,"Right scaling vector wrong length"); 1299 VecGetArray_Fast(rr,r); 1300 v = a->a; jj = a->j; 1301 for ( i=0; i<nz; i++ ) { 1302 (*v++) *= r[*jj++ + shift]; 1303 } 1304 PLogFlops(nz); 1305 } 1306 return 0; 1307 } 1308 1309 #undef __FUNC__ 1310 #define __FUNC__ "MatGetSubMatrix_SeqAIJ" 1311 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,MatGetSubMatrixCall scall,Mat *B) 1312 { 1313 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1314 int nznew, *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1315 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1316 register int sum,lensi; 1317 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1318 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1319 Scalar *a_new,*mat_a; 1320 Mat C; 1321 1322 ierr = ISSorted(isrow,(PetscTruth*)&i); 1323 if (!i) SETERRQ(1,0,"ISrow is not sorted"); 1324 ierr = ISSorted(iscol,(PetscTruth*)&i); 1325 if (!i) SETERRQ(1,0,"IScol is not sorted"); 1326 1327 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1328 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1329 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1330 1331 if (ISStrideGetInfo(iscol,&first,&step) && step == 1) { /* no need to sort */ 1332 /* special case of contiguous rows */ 1333 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1334 starts = lens + ncols; 1335 /* loop over new rows determining lens and starting points */ 1336 for (i=0; i<nrows; i++) { 1337 kstart = ai[irow[i]]+shift; 1338 kend = kstart + ailen[irow[i]]; 1339 for ( k=kstart; k<kend; k++ ) { 1340 if (aj[k]+shift >= first) { 1341 starts[i] = k; 1342 break; 1343 } 1344 } 1345 sum = 0; 1346 while (k < kend) { 1347 if (aj[k++]+shift >= first+ncols) break; 1348 sum++; 1349 } 1350 lens[i] = sum; 1351 } 1352 /* create submatrix */ 1353 if (scall == MAT_REUSE_MATRIX) { 1354 int n_cols,n_rows; 1355 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1356 if (n_rows != nrows || n_cols != ncols) SETERRQ(1,0,""); 1357 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1358 C = *B; 1359 } 1360 else { 1361 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1362 } 1363 c = (Mat_SeqAIJ*) C->data; 1364 1365 /* loop over rows inserting into submatrix */ 1366 a_new = c->a; 1367 j_new = c->j; 1368 i_new = c->i; 1369 i_new[0] = -shift; 1370 for (i=0; i<nrows; i++) { 1371 ii = starts[i]; 1372 lensi = lens[i]; 1373 for ( k=0; k<lensi; k++ ) { 1374 *j_new++ = aj[ii+k] - first; 1375 } 1376 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1377 a_new += lensi; 1378 i_new[i+1] = i_new[i] + lensi; 1379 c->ilen[i] = lensi; 1380 } 1381 PetscFree(lens); 1382 } 1383 else { 1384 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1385 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1386 ssmap = smap + shift; 1387 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1388 PetscMemzero(smap,oldcols*sizeof(int)); 1389 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1390 /* determine lens of each row */ 1391 for (i=0; i<nrows; i++) { 1392 kstart = ai[irow[i]]+shift; 1393 kend = kstart + a->ilen[irow[i]]; 1394 lens[i] = 0; 1395 for ( k=kstart; k<kend; k++ ) { 1396 if (ssmap[aj[k]]) { 1397 lens[i]++; 1398 } 1399 } 1400 } 1401 /* Create and fill new matrix */ 1402 if (scall == MAT_REUSE_MATRIX) { 1403 c = (Mat_SeqAIJ *)((*B)->data); 1404 1405 if (c->m != nrows || c->n != ncols) SETERRQ(1,0,""); 1406 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1407 SETERRQ(1,0,"Cannot reuse matrix. wrong no of nonzeros"); 1408 } 1409 PetscMemzero(c->ilen,c->m*sizeof(int)); 1410 C = *B; 1411 } 1412 else { 1413 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1414 } 1415 c = (Mat_SeqAIJ *)(C->data); 1416 for (i=0; i<nrows; i++) { 1417 row = irow[i]; 1418 nznew = 0; 1419 kstart = ai[row]+shift; 1420 kend = kstart + a->ilen[row]; 1421 mat_i = c->i[i]+shift; 1422 mat_j = c->j + mat_i; 1423 mat_a = c->a + mat_i; 1424 mat_ilen = c->ilen + i; 1425 for ( k=kstart; k<kend; k++ ) { 1426 if ((tcol=ssmap[a->j[k]])) { 1427 *mat_j++ = tcol - (!shift); 1428 *mat_a++ = a->a[k]; 1429 (*mat_ilen)++; 1430 1431 } 1432 } 1433 } 1434 /* Free work space */ 1435 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1436 PetscFree(smap); PetscFree(lens); 1437 } 1438 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1439 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1440 1441 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1442 *B = C; 1443 return 0; 1444 } 1445 1446 /* 1447 note: This can only work for identity for row and col. It would 1448 be good to check this and otherwise generate an error. 1449 */ 1450 #undef __FUNC__ 1451 #define __FUNC__ "MatILUFactor_SeqAIJ" 1452 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1453 { 1454 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1455 int ierr; 1456 Mat outA; 1457 1458 if (fill != 0) SETERRQ(1,0,"Only fill=0 supported"); 1459 1460 outA = inA; 1461 inA->factor = FACTOR_LU; 1462 a->row = row; 1463 a->col = col; 1464 1465 if (!a->solve_work) { /* this matrix may have been factored before */ 1466 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar)); CHKPTRQ(a->solve_work); 1467 } 1468 1469 if (!a->diag) { 1470 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1471 } 1472 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1473 return 0; 1474 } 1475 1476 #include "pinclude/plapack.h" 1477 #undef __FUNC__ 1478 #define __FUNC__ "MatScale_SeqAIJ" 1479 int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1480 { 1481 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1482 int one = 1; 1483 BLscal_( &a->nz, alpha, a->a, &one ); 1484 PLogFlops(a->nz); 1485 return 0; 1486 } 1487 1488 #undef __FUNC__ 1489 #define __FUNC__ "MatGetSubMatrices_SeqAIJ" 1490 int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall, 1491 Mat **B) 1492 { 1493 int ierr,i; 1494 1495 if (scall == MAT_INITIAL_MATRIX) { 1496 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1497 } 1498 1499 for ( i=0; i<n; i++ ) { 1500 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 1501 } 1502 return 0; 1503 } 1504 1505 #undef __FUNC__ 1506 #define __FUNC__ "MatGetBlockSize_SeqAIJ" 1507 int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1508 { 1509 *bs = 1; 1510 return 0; 1511 } 1512 1513 #undef __FUNC__ 1514 #define __FUNC__ "MatIncreaseOverlap_SeqAIJ" 1515 int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1516 { 1517 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1518 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1519 int start, end, *ai, *aj; 1520 char *table; 1521 shift = a->indexshift; 1522 m = a->m; 1523 ai = a->i; 1524 aj = a->j+shift; 1525 1526 if (ov < 0) SETERRQ(1,0,"illegal overlap value used"); 1527 1528 table = (char *) PetscMalloc((m/BITSPERBYTE +1)*sizeof(char)); CHKPTRQ(table); 1529 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1530 1531 for ( i=0; i<is_max; i++ ) { 1532 /* Initialize the two local arrays */ 1533 isz = 0; 1534 PetscMemzero(table,(m/BITSPERBYTE +1)*sizeof(char)); 1535 1536 /* Extract the indices, assume there can be duplicate entries */ 1537 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1538 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1539 1540 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1541 for ( j=0; j<n ; ++j){ 1542 if(!BT_LOOKUP(table, idx[j])) { nidx[isz++] = idx[j];} 1543 } 1544 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1545 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1546 1547 k = 0; 1548 for ( j=0; j<ov; j++){ /* for each overlap */ 1549 n = isz; 1550 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1551 row = nidx[k]; 1552 start = ai[row]; 1553 end = ai[row+1]; 1554 for ( l = start; l<end ; l++){ 1555 val = aj[l] + shift; 1556 if (!BT_LOOKUP(table,val)) {nidx[isz++] = val;} 1557 } 1558 } 1559 } 1560 ierr = ISCreateGeneral(PETSC_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1561 } 1562 PetscFree(table); 1563 PetscFree(nidx); 1564 return 0; 1565 } 1566 1567 /* -------------------------------------------------------------- */ 1568 #undef __FUNC__ 1569 #define __FUNC__ "MatPermute_SeqAIJ" 1570 int MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 1571 { 1572 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1573 Scalar *vwork; 1574 int i, ierr, nz, m = a->m, n = a->n, *cwork; 1575 int *row,*col,*cnew,j,*lens; 1576 IS icolp,irowp; 1577 1578 ierr = ISInvertPermutation(rowp,&irowp); CHKERRQ(ierr); 1579 ierr = ISGetIndices(irowp,&row); CHKERRQ(ierr); 1580 ierr = ISInvertPermutation(colp,&icolp); CHKERRQ(ierr); 1581 ierr = ISGetIndices(icolp,&col); CHKERRQ(ierr); 1582 1583 /* determine lengths of permuted rows */ 1584 lens = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(lens); 1585 for (i=0; i<m; i++ ) { 1586 lens[row[i]] = a->i[i+1] - a->i[i]; 1587 } 1588 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1589 PetscFree(lens); 1590 1591 cnew = (int *) PetscMalloc( n*sizeof(int) ); CHKPTRQ(cnew); 1592 for (i=0; i<m; i++) { 1593 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1594 for (j=0; j<nz; j++ ) { cnew[j] = col[cwork[j]];} 1595 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES); CHKERRQ(ierr); 1596 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1597 } 1598 PetscFree(cnew); 1599 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1600 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1601 ierr = ISRestoreIndices(irowp,&row); CHKERRQ(ierr); 1602 ierr = ISRestoreIndices(icolp,&col); CHKERRQ(ierr); 1603 ierr = ISDestroy(irowp); CHKERRQ(ierr); 1604 ierr = ISDestroy(icolp); CHKERRQ(ierr); 1605 return 0; 1606 } 1607 1608 #undef __FUNC__ 1609 #define __FUNC__ "MatPrintHelp_SeqAIJ" 1610 int MatPrintHelp_SeqAIJ(Mat A) 1611 { 1612 static int called = 0; 1613 MPI_Comm comm = A->comm; 1614 1615 if (called) return 0; else called = 1; 1616 PetscPrintf(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1617 PetscPrintf(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n"); 1618 PetscPrintf(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n"); 1619 PetscPrintf(comm," -mat_aij_no_inode: Do not use inodes\n"); 1620 PetscPrintf(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n"); 1621 #if defined(HAVE_ESSL) 1622 PetscPrintf(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n"); 1623 #endif 1624 return 0; 1625 } 1626 extern int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1627 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1628 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1629 1630 /* -------------------------------------------------------------------*/ 1631 static struct _MatOps MatOps = {MatSetValues_SeqAIJ, 1632 MatGetRow_SeqAIJ,MatRestoreRow_SeqAIJ, 1633 MatMult_SeqAIJ,MatMultAdd_SeqAIJ, 1634 MatMultTrans_SeqAIJ,MatMultTransAdd_SeqAIJ, 1635 MatSolve_SeqAIJ,MatSolveAdd_SeqAIJ, 1636 MatSolveTrans_SeqAIJ,MatSolveTransAdd_SeqAIJ, 1637 MatLUFactor_SeqAIJ,0, 1638 MatRelax_SeqAIJ, 1639 MatTranspose_SeqAIJ, 1640 MatGetInfo_SeqAIJ,MatEqual_SeqAIJ, 1641 MatGetDiagonal_SeqAIJ,MatDiagonalScale_SeqAIJ,MatNorm_SeqAIJ, 1642 0,MatAssemblyEnd_SeqAIJ, 1643 MatCompress_SeqAIJ, 1644 MatSetOption_SeqAIJ,MatZeroEntries_SeqAIJ,MatZeroRows_SeqAIJ, 1645 MatLUFactorSymbolic_SeqAIJ,MatLUFactorNumeric_SeqAIJ,0,0, 1646 MatGetSize_SeqAIJ,MatGetSize_SeqAIJ,MatGetOwnershipRange_SeqAIJ, 1647 MatILUFactorSymbolic_SeqAIJ,0, 1648 0,0, 1649 MatConvertSameType_SeqAIJ,0,0, 1650 MatILUFactor_SeqAIJ,0,0, 1651 MatGetSubMatrices_SeqAIJ,MatIncreaseOverlap_SeqAIJ, 1652 MatGetValues_SeqAIJ,0, 1653 MatPrintHelp_SeqAIJ, 1654 MatScale_SeqAIJ,0,0, 1655 MatILUDTFactor_SeqAIJ, 1656 MatGetBlockSize_SeqAIJ, 1657 MatGetRowIJ_SeqAIJ, 1658 MatRestoreRowIJ_SeqAIJ, 1659 MatGetColumnIJ_SeqAIJ, 1660 MatRestoreColumnIJ_SeqAIJ, 1661 MatFDColoringCreate_SeqAIJ, 1662 MatColoringPatch_SeqAIJ, 1663 0, 1664 MatPermute_SeqAIJ}; 1665 1666 extern int MatUseSuperLU_SeqAIJ(Mat); 1667 extern int MatUseEssl_SeqAIJ(Mat); 1668 extern int MatUseDXML_SeqAIJ(Mat); 1669 1670 #undef __FUNC__ 1671 #define __FUNC__ "MatCreateSeqAIJ" 1672 /*@C 1673 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 1674 (the default parallel PETSc format). For good matrix assembly performance 1675 the user should preallocate the matrix storage by setting the parameter nz 1676 (or the array nzz). By setting these parameters accurately, performance 1677 during matrix assembly can be increased by more than a factor of 50. 1678 1679 Input Parameters: 1680 . comm - MPI communicator, set to PETSC_COMM_SELF 1681 . m - number of rows 1682 . n - number of columns 1683 . nz - number of nonzeros per row (same for all rows) 1684 . nzz - array containing the number of nonzeros in the various rows 1685 (possibly different for each row) or PETSC_NULL 1686 1687 Output Parameter: 1688 . A - the matrix 1689 1690 Notes: 1691 The AIJ format (also called the Yale sparse matrix format or 1692 compressed row storage), is fully compatible with standard Fortran 77 1693 storage. That is, the stored row and column indices can begin at 1694 either one (as in Fortran) or zero. See the users' manual for details. 1695 1696 Specify the preallocated storage with either nz or nnz (not both). 1697 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1698 allocation. For large problems you MUST preallocate memory or you 1699 will get TERRIBLE performance, see the users' manual chapter on matrices. 1700 1701 By default, this format uses inodes (identical nodes) when possible, to 1702 improve numerical efficiency of matrix-vector products and solves. We 1703 search for consecutive rows with the same nonzero structure, thereby 1704 reusing matrix information to achieve increased efficiency. 1705 1706 Options Database Keys: 1707 $ -mat_aij_no_inode - Do not use inodes 1708 $ -mat_aij_inode_limit <limit> - Set inode limit. 1709 $ (max limit=5) 1710 $ -mat_aij_oneindex - Internally use indexing starting at 1 1711 $ rather than 0. Note: When calling MatSetValues(), 1712 $ the user still MUST index entries starting at 0! 1713 1714 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues() 1715 @*/ 1716 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 1717 { 1718 Mat B; 1719 Mat_SeqAIJ *b; 1720 int i, len, ierr, flg,size; 1721 1722 MPI_Comm_size(comm,&size); 1723 if (size > 1) SETERRQ(1,0,"Comm must be of size 1"); 1724 1725 *A = 0; 1726 PetscHeaderCreate(B,_p_Mat,MAT_COOKIE,MATSEQAIJ,comm,MatDestroy,MatView); 1727 PLogObjectCreate(B); 1728 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 1729 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 1730 PetscMemcpy(&B->ops,&MatOps,sizeof(struct _MatOps)); 1731 B->destroy = MatDestroy_SeqAIJ; 1732 B->view = MatView_SeqAIJ; 1733 B->factor = 0; 1734 B->lupivotthreshold = 1.0; 1735 B->mapping = 0; 1736 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold, 1737 &flg); CHKERRQ(ierr); 1738 b->ilu_preserve_row_sums = PETSC_FALSE; 1739 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums", 1740 (int*) &b->ilu_preserve_row_sums); CHKERRQ(ierr); 1741 b->row = 0; 1742 b->col = 0; 1743 b->indexshift = 0; 1744 b->reallocs = 0; 1745 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 1746 if (flg) b->indexshift = -1; 1747 1748 b->m = m; B->m = m; B->M = m; 1749 b->n = n; B->n = n; B->N = n; 1750 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 1751 if (nnz == PETSC_NULL) { 1752 if (nz == PETSC_DEFAULT) nz = 10; 1753 else if (nz <= 0) nz = 1; 1754 for ( i=0; i<m; i++ ) b->imax[i] = nz; 1755 nz = nz*m; 1756 } 1757 else { 1758 nz = 0; 1759 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1760 } 1761 1762 /* allocate the matrix space */ 1763 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 1764 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1765 b->j = (int *) (b->a + nz); 1766 PetscMemzero(b->j,nz*sizeof(int)); 1767 b->i = b->j + nz; 1768 b->singlemalloc = 1; 1769 1770 b->i[0] = -b->indexshift; 1771 for (i=1; i<m+1; i++) { 1772 b->i[i] = b->i[i-1] + b->imax[i-1]; 1773 } 1774 1775 /* b->ilen will count nonzeros in each row so far. */ 1776 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 1777 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1778 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 1779 1780 b->nz = 0; 1781 b->maxnz = nz; 1782 b->sorted = 0; 1783 b->roworiented = 1; 1784 b->nonew = 0; 1785 b->diag = 0; 1786 b->solve_work = 0; 1787 b->spptr = 0; 1788 b->inode.node_count = 0; 1789 b->inode.size = 0; 1790 b->inode.limit = 5; 1791 b->inode.max_limit = 5; 1792 B->info.nz_unneeded = (double)b->maxnz; 1793 1794 *A = B; 1795 1796 /* SuperLU is not currently supported through PETSc */ 1797 #if defined(HAVE_SUPERLU) 1798 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 1799 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 1800 #endif 1801 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 1802 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 1803 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 1804 if (flg) { 1805 if (!b->indexshift) SETERRQ(1,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 1806 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 1807 } 1808 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1809 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1810 return 0; 1811 } 1812 1813 #undef __FUNC__ 1814 #define __FUNC__ "MatConvertSameType_SeqAIJ" 1815 int MatConvertSameType_SeqAIJ(Mat A,Mat *B,int cpvalues) 1816 { 1817 Mat C; 1818 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 1819 int i,len, m = a->m,shift = a->indexshift; 1820 1821 *B = 0; 1822 PetscHeaderCreate(C,_p_Mat,MAT_COOKIE,MATSEQAIJ,A->comm,MatDestroy,MatView); 1823 PLogObjectCreate(C); 1824 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 1825 PetscMemcpy(&C->ops,&A->ops,sizeof(struct _MatOps)); 1826 C->destroy = MatDestroy_SeqAIJ; 1827 C->view = MatView_SeqAIJ; 1828 C->factor = A->factor; 1829 c->row = 0; 1830 c->col = 0; 1831 c->indexshift = shift; 1832 C->assembled = PETSC_TRUE; 1833 1834 c->m = C->m = a->m; 1835 c->n = C->n = a->n; 1836 C->M = a->m; 1837 C->N = a->n; 1838 1839 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 1840 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 1841 for ( i=0; i<m; i++ ) { 1842 c->imax[i] = a->imax[i]; 1843 c->ilen[i] = a->ilen[i]; 1844 } 1845 1846 /* allocate the matrix space */ 1847 c->singlemalloc = 1; 1848 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 1849 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1850 c->j = (int *) (c->a + a->i[m] + shift); 1851 c->i = c->j + a->i[m] + shift; 1852 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 1853 if (m > 0) { 1854 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 1855 if (cpvalues == COPY_VALUES) { 1856 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 1857 } 1858 } 1859 1860 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1861 c->sorted = a->sorted; 1862 c->roworiented = a->roworiented; 1863 c->nonew = a->nonew; 1864 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 1865 1866 if (a->diag) { 1867 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 1868 PLogObjectMemory(C,(m+1)*sizeof(int)); 1869 for ( i=0; i<m; i++ ) { 1870 c->diag[i] = a->diag[i]; 1871 } 1872 } 1873 else c->diag = 0; 1874 c->inode.limit = a->inode.limit; 1875 c->inode.max_limit = a->inode.max_limit; 1876 if (a->inode.size){ 1877 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 1878 c->inode.node_count = a->inode.node_count; 1879 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 1880 } else { 1881 c->inode.size = 0; 1882 c->inode.node_count = 0; 1883 } 1884 c->nz = a->nz; 1885 c->maxnz = a->maxnz; 1886 c->solve_work = 0; 1887 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1888 1889 *B = C; 1890 return 0; 1891 } 1892 1893 #undef __FUNC__ 1894 #define __FUNC__ "MatLoad_SeqAIJ" 1895 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 1896 { 1897 Mat_SeqAIJ *a; 1898 Mat B; 1899 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 1900 MPI_Comm comm; 1901 1902 PetscObjectGetComm((PetscObject) viewer,&comm); 1903 MPI_Comm_size(comm,&size); 1904 if (size > 1) SETERRQ(1,0,"view must have one processor"); 1905 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 1906 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 1907 if (header[0] != MAT_COOKIE) SETERRQ(1,0,"not matrix object in file"); 1908 M = header[1]; N = header[2]; nz = header[3]; 1909 1910 /* read in row lengths */ 1911 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 1912 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 1913 1914 /* create our matrix */ 1915 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 1916 B = *A; 1917 a = (Mat_SeqAIJ *) B->data; 1918 shift = a->indexshift; 1919 1920 /* read in column indices and adjust for Fortran indexing*/ 1921 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT); CHKERRQ(ierr); 1922 if (shift) { 1923 for ( i=0; i<nz; i++ ) { 1924 a->j[i] += 1; 1925 } 1926 } 1927 1928 /* read in nonzero values */ 1929 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR); CHKERRQ(ierr); 1930 1931 /* set matrix "i" values */ 1932 a->i[0] = -shift; 1933 for ( i=1; i<= M; i++ ) { 1934 a->i[i] = a->i[i-1] + rowlengths[i-1]; 1935 a->ilen[i-1] = rowlengths[i-1]; 1936 } 1937 PetscFree(rowlengths); 1938 1939 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1940 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1941 return 0; 1942 } 1943 1944 #undef __FUNC__ 1945 #define __FUNC__ "MatEqual_SeqAIJ" 1946 int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 1947 { 1948 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 1949 1950 if (B->type !=MATSEQAIJ)SETERRQ(1,0,"Matrices must be same type"); 1951 1952 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 1953 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 1954 (a->indexshift != b->indexshift)) { 1955 *flg = PETSC_FALSE; return 0; 1956 } 1957 1958 /* if the a->i are the same */ 1959 if (PetscMemcmp(a->i,b->i,(a->m+1)*sizeof(int))) { 1960 *flg = PETSC_FALSE; return 0; 1961 } 1962 1963 /* if a->j are the same */ 1964 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 1965 *flg = PETSC_FALSE; return 0; 1966 } 1967 1968 /* if a->a are the same */ 1969 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 1970 *flg = PETSC_FALSE; return 0; 1971 } 1972 *flg = PETSC_TRUE; 1973 return 0; 1974 1975 } 1976