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