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