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