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