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