1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: aij.c,v 1.264 1998/05/13 14:13:43 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 A->qlist = 0; 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/blaslapack.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__ "MatSeqAIJSetColumnIndices_SeqAIJ" 1755 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 1756 { 1757 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1758 int i,nz,n; 1759 1760 PetscFunctionBegin; 1761 if (aij->indexshift) SETERRQ(1,1,"No support with 1 based indexing"); 1762 1763 nz = aij->maxnz; 1764 n = aij->n; 1765 for (i=0; i<nz; i++) { 1766 aij->j[i] = indices[i]; 1767 } 1768 aij->nz = nz; 1769 for ( i=0; i<n; i++ ) { 1770 aij->ilen[i] = aij->imax[i]; 1771 } 1772 1773 PetscFunctionReturn(0); 1774 } 1775 1776 #undef __FUNC__ 1777 #define __FUNC__ "MatSeqAIJSetColumnIndices" 1778 /*@ 1779 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 1780 in the matrix. 1781 1782 Input Parameters: 1783 + mat - the SeqAIJ matrix 1784 - indices - the column indices 1785 1786 Notes: 1787 This can be called if you have precomputed the nonzero structure of the 1788 matrix and want to provide it to the matrix object to improve the performance 1789 of the MatSetValues() operation. 1790 1791 You MUST have set the correct numbers of nonzeros per row in the call to 1792 MatCreateSeqAIJ(). 1793 1794 MUST be called before any calls to MatSetValues(); 1795 1796 @*/ 1797 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 1798 { 1799 int ierr,(*f)(Mat,int *); 1800 1801 PetscFunctionBegin; 1802 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1803 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void **)&f); 1804 CHKERRQ(ierr); 1805 if (f) { 1806 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1807 } else { 1808 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1809 } 1810 PetscFunctionReturn(0); 1811 } 1812 1813 #undef __FUNC__ 1814 #define __FUNC__ "MatCreateSeqAIJ" 1815 /*@C 1816 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 1817 (the default parallel PETSc format). For good matrix assembly performance 1818 the user should preallocate the matrix storage by setting the parameter nz 1819 (or the array nzz). By setting these parameters accurately, performance 1820 during matrix assembly can be increased by more than a factor of 50. 1821 1822 Collective on MPI_Comm 1823 1824 Input Parameters: 1825 + comm - MPI communicator, set to PETSC_COMM_SELF 1826 . m - number of rows 1827 . n - number of columns 1828 . nz - number of nonzeros per row (same for all rows) 1829 - nzz - array containing the number of nonzeros in the various rows 1830 (possibly different for each row) or PETSC_NULL 1831 1832 Output Parameter: 1833 . A - the matrix 1834 1835 Notes: 1836 The AIJ format (also called the Yale sparse matrix format or 1837 compressed row storage), is fully compatible with standard Fortran 77 1838 storage. That is, the stored row and column indices can begin at 1839 either one (as in Fortran) or zero. See the users' manual for details. 1840 1841 Specify the preallocated storage with either nz or nnz (not both). 1842 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1843 allocation. For large problems you MUST preallocate memory or you 1844 will get TERRIBLE performance, see the users' manual chapter on matrices. 1845 1846 By default, this format uses inodes (identical nodes) when possible, to 1847 improve numerical efficiency of matrix-vector products and solves. We 1848 search for consecutive rows with the same nonzero structure, thereby 1849 reusing matrix information to achieve increased efficiency. 1850 1851 Options Database Keys: 1852 + -mat_aij_no_inode - Do not use inodes 1853 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 1854 - -mat_aij_oneindex - Internally use indexing starting at 1 1855 rather than 0. Note that when calling MatSetValues(), 1856 the user still MUST index entries starting at 0! 1857 1858 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices() 1859 @*/ 1860 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 1861 { 1862 Mat B; 1863 Mat_SeqAIJ *b; 1864 int i, len, ierr, flg,size; 1865 1866 PetscFunctionBegin; 1867 MPI_Comm_size(comm,&size); 1868 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Comm must be of size 1"); 1869 1870 *A = 0; 1871 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,comm,MatDestroy,MatView); 1872 PLogObjectCreate(B); 1873 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 1874 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 1875 PetscMemcpy(B->ops,&MatOps,sizeof(struct _MatOps)); 1876 B->ops->destroy = MatDestroy_SeqAIJ; 1877 B->ops->view = MatView_SeqAIJ; 1878 B->factor = 0; 1879 B->lupivotthreshold = 1.0; 1880 B->mapping = 0; 1881 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,&flg); CHKERRQ(ierr); 1882 b->ilu_preserve_row_sums = PETSC_FALSE; 1883 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",(int*) &b->ilu_preserve_row_sums);CHKERRQ(ierr); 1884 b->row = 0; 1885 b->col = 0; 1886 b->icol = 0; 1887 b->indexshift = 0; 1888 b->reallocs = 0; 1889 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 1890 if (flg) b->indexshift = -1; 1891 1892 b->m = m; B->m = m; B->M = m; 1893 b->n = n; B->n = n; B->N = n; 1894 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 1895 if (nnz == PETSC_NULL) { 1896 if (nz == PETSC_DEFAULT) nz = 10; 1897 else if (nz <= 0) nz = 1; 1898 for ( i=0; i<m; i++ ) b->imax[i] = nz; 1899 nz = nz*m; 1900 } else { 1901 nz = 0; 1902 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1903 } 1904 1905 /* allocate the matrix space */ 1906 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 1907 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1908 b->j = (int *) (b->a + nz); 1909 PetscMemzero(b->j,nz*sizeof(int)); 1910 b->i = b->j + nz; 1911 b->singlemalloc = 1; 1912 1913 b->i[0] = -b->indexshift; 1914 for (i=1; i<m+1; i++) { 1915 b->i[i] = b->i[i-1] + b->imax[i-1]; 1916 } 1917 1918 /* b->ilen will count nonzeros in each row so far. */ 1919 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 1920 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1921 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 1922 1923 b->nz = 0; 1924 b->maxnz = nz; 1925 b->sorted = 0; 1926 b->roworiented = 1; 1927 b->nonew = 0; 1928 b->diag = 0; 1929 b->solve_work = 0; 1930 b->spptr = 0; 1931 b->inode.node_count = 0; 1932 b->inode.size = 0; 1933 b->inode.limit = 5; 1934 b->inode.max_limit = 5; 1935 B->info.nz_unneeded = (double)b->maxnz; 1936 1937 *A = B; 1938 1939 /* SuperLU is not currently supported through PETSc */ 1940 #if defined(HAVE_SUPERLU) 1941 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 1942 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 1943 #endif 1944 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 1945 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 1946 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 1947 if (flg) { 1948 if (!b->indexshift) SETERRQ( PETSC_ERR_LIB,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 1949 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 1950 } 1951 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1952 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1953 1954 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 1955 "MatSeqAIJSetColumnIndices_SeqAIJ", 1956 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 1957 PetscFunctionReturn(0); 1958 } 1959 1960 #undef __FUNC__ 1961 #define __FUNC__ "MatConvertSameType_SeqAIJ" 1962 int MatConvertSameType_SeqAIJ(Mat A,Mat *B,int cpvalues) 1963 { 1964 Mat C; 1965 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 1966 int i,len, m = a->m,shift = a->indexshift,ierr; 1967 1968 PetscFunctionBegin; 1969 *B = 0; 1970 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,A->comm,MatDestroy,MatView); 1971 PLogObjectCreate(C); 1972 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 1973 PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps)); 1974 C->ops->destroy = MatDestroy_SeqAIJ; 1975 C->ops->view = MatView_SeqAIJ; 1976 C->factor = A->factor; 1977 c->row = 0; 1978 c->col = 0; 1979 c->icol = 0; 1980 c->indexshift = shift; 1981 C->assembled = PETSC_TRUE; 1982 1983 c->m = C->m = a->m; 1984 c->n = C->n = a->n; 1985 C->M = a->m; 1986 C->N = a->n; 1987 1988 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 1989 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 1990 for ( i=0; i<m; i++ ) { 1991 c->imax[i] = a->imax[i]; 1992 c->ilen[i] = a->ilen[i]; 1993 } 1994 1995 /* allocate the matrix space */ 1996 c->singlemalloc = 1; 1997 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 1998 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1999 c->j = (int *) (c->a + a->i[m] + shift); 2000 c->i = c->j + a->i[m] + shift; 2001 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 2002 if (m > 0) { 2003 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 2004 if (cpvalues == COPY_VALUES) { 2005 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 2006 } 2007 } 2008 2009 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2010 c->sorted = a->sorted; 2011 c->roworiented = a->roworiented; 2012 c->nonew = a->nonew; 2013 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2014 2015 if (a->diag) { 2016 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 2017 PLogObjectMemory(C,(m+1)*sizeof(int)); 2018 for ( i=0; i<m; i++ ) { 2019 c->diag[i] = a->diag[i]; 2020 } 2021 } else c->diag = 0; 2022 c->inode.limit = a->inode.limit; 2023 c->inode.max_limit = a->inode.max_limit; 2024 if (a->inode.size){ 2025 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 2026 c->inode.node_count = a->inode.node_count; 2027 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 2028 } else { 2029 c->inode.size = 0; 2030 c->inode.node_count = 0; 2031 } 2032 c->nz = a->nz; 2033 c->maxnz = a->maxnz; 2034 c->solve_work = 0; 2035 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2036 2037 *B = C; 2038 ierr = PetscObjectComposeFunction((PetscObject)C,"MatSeqAIJSetColumnIndices_C", 2039 "MatSeqAIJSetColumnIndices_SeqAIJ", 2040 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2041 PetscFunctionReturn(0); 2042 } 2043 2044 #undef __FUNC__ 2045 #define __FUNC__ "MatLoad_SeqAIJ" 2046 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 2047 { 2048 Mat_SeqAIJ *a; 2049 Mat B; 2050 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 2051 MPI_Comm comm; 2052 2053 PetscFunctionBegin; 2054 PetscObjectGetComm((PetscObject) viewer,&comm); 2055 MPI_Comm_size(comm,&size); 2056 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,0,"view must have one processor"); 2057 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 2058 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 2059 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not matrix object in file"); 2060 M = header[1]; N = header[2]; nz = header[3]; 2061 2062 if (nz < 0) { 2063 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 2064 } 2065 2066 /* read in row lengths */ 2067 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 2068 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 2069 2070 /* create our matrix */ 2071 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 2072 B = *A; 2073 a = (Mat_SeqAIJ *) B->data; 2074 shift = a->indexshift; 2075 2076 /* read in column indices and adjust for Fortran indexing*/ 2077 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT); CHKERRQ(ierr); 2078 if (shift) { 2079 for ( i=0; i<nz; i++ ) { 2080 a->j[i] += 1; 2081 } 2082 } 2083 2084 /* read in nonzero values */ 2085 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR); CHKERRQ(ierr); 2086 2087 /* set matrix "i" values */ 2088 a->i[0] = -shift; 2089 for ( i=1; i<= M; i++ ) { 2090 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2091 a->ilen[i-1] = rowlengths[i-1]; 2092 } 2093 PetscFree(rowlengths); 2094 2095 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2096 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2097 PetscFunctionReturn(0); 2098 } 2099 2100 #undef __FUNC__ 2101 #define __FUNC__ "MatEqual_SeqAIJ" 2102 int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 2103 { 2104 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 2105 2106 PetscFunctionBegin; 2107 if (B->type !=MATSEQAIJ)SETERRQ(PETSC_ERR_ARG_INCOMP,0,"Matrices must be same type"); 2108 2109 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 2110 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 2111 (a->indexshift != b->indexshift)) { 2112 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2113 } 2114 2115 /* if the a->i are the same */ 2116 if (PetscMemcmp(a->i,b->i,(a->m+1)*sizeof(int))) { 2117 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2118 } 2119 2120 /* if a->j are the same */ 2121 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 2122 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2123 } 2124 2125 /* if a->a are the same */ 2126 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 2127 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2128 } 2129 *flg = PETSC_TRUE; 2130 PetscFunctionReturn(0); 2131 2132 } 2133