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