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