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