1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: aij.c,v 1.280 1998/08/03 14:58:58 balay 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 "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 diag[i] = a->i[i+1]; 916 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 917 if (a->j[j]+shift == i) { 918 diag[i] = j - shift; 919 break; 920 } 921 } 922 } 923 a->diag = diag; 924 PetscFunctionReturn(0); 925 } 926 927 #undef __FUNC__ 928 #define __FUNC__ "MatRelax_SeqAIJ" 929 int MatRelax_SeqAIJ(Mat A,Vec bb,double omega,MatSORType flag, 930 double fshift,int its,Vec xx) 931 { 932 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 933 Scalar *x, *b, *bs, d, *xs, sum, *v = a->a,*t,scale,*ts, *xb; 934 int ierr, *idx, *diag,n = a->n, m = a->m, i, shift = a->indexshift; 935 936 PetscFunctionBegin; 937 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 938 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 939 if (!a->diag) {ierr = MatMarkDiag_SeqAIJ(A);CHKERRQ(ierr);} 940 diag = a->diag; 941 xs = x + shift; /* shifted by one for index start of a or a->j*/ 942 if (flag == SOR_APPLY_UPPER) { 943 /* apply ( U + D/omega) to the vector */ 944 bs = b + shift; 945 for ( i=0; i<m; i++ ) { 946 d = fshift + a->a[diag[i] + shift]; 947 n = a->i[i+1] - diag[i] - 1; 948 PLogFlops(2*n-1); 949 idx = a->j + diag[i] + (!shift); 950 v = a->a + diag[i] + (!shift); 951 sum = b[i]*d/omega; 952 SPARSEDENSEDOT(sum,bs,v,idx,n); 953 x[i] = sum; 954 } 955 PetscFunctionReturn(0); 956 } 957 if (flag == SOR_APPLY_LOWER) { 958 SETERRQ(PETSC_ERR_SUP,0,"SOR_APPLY_LOWER is not done"); 959 } else if (flag & SOR_EISENSTAT) { 960 /* Let A = L + U + D; where L is lower trianglar, 961 U is upper triangular, E is diagonal; This routine applies 962 963 (L + E)^{-1} A (U + E)^{-1} 964 965 to a vector efficiently using Eisenstat's trick. This is for 966 the case of SSOR preconditioner, so E is D/omega where omega 967 is the relaxation factor. 968 */ 969 t = (Scalar *) PetscMalloc( m*sizeof(Scalar) ); CHKPTRQ(t); 970 scale = (2.0/omega) - 1.0; 971 972 /* x = (E + U)^{-1} b */ 973 for ( i=m-1; i>=0; i-- ) { 974 d = fshift + a->a[diag[i] + shift]; 975 n = a->i[i+1] - diag[i] - 1; 976 PLogFlops(2*n-1); 977 idx = a->j + diag[i] + (!shift); 978 v = a->a + diag[i] + (!shift); 979 sum = b[i]; 980 SPARSEDENSEMDOT(sum,xs,v,idx,n); 981 x[i] = omega*(sum/d); 982 } 983 984 /* t = b - (2*E - D)x */ 985 v = a->a; 986 for ( i=0; i<m; i++ ) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 987 PLogFlops(2*m); 988 989 990 /* t = (E + L)^{-1}t */ 991 ts = t + shift; /* shifted by one for index start of a or a->j*/ 992 diag = a->diag; 993 for ( i=0; i<m; i++ ) { 994 d = fshift + a->a[diag[i]+shift]; 995 n = diag[i] - a->i[i]; 996 PLogFlops(2*n-1); 997 idx = a->j + a->i[i] + shift; 998 v = a->a + a->i[i] + shift; 999 sum = t[i]; 1000 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1001 t[i] = omega*(sum/d); 1002 } 1003 1004 /* x = x + t */ 1005 for ( i=0; i<m; i++ ) { x[i] += t[i]; } 1006 PetscFree(t); 1007 PetscFunctionReturn(0); 1008 } 1009 if (flag & SOR_ZERO_INITIAL_GUESS) { 1010 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1011 for ( i=0; i<m; i++ ) { 1012 d = fshift + a->a[diag[i]+shift]; 1013 n = diag[i] - a->i[i]; 1014 PLogFlops(2*n-1); 1015 idx = a->j + a->i[i] + shift; 1016 v = a->a + a->i[i] + shift; 1017 sum = b[i]; 1018 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1019 x[i] = omega*(sum/d); 1020 } 1021 xb = x; 1022 } else xb = b; 1023 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1024 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1025 for ( i=0; i<m; i++ ) { 1026 x[i] *= a->a[diag[i]+shift]; 1027 } 1028 PLogFlops(m); 1029 } 1030 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1031 for ( i=m-1; i>=0; i-- ) { 1032 d = fshift + a->a[diag[i] + shift]; 1033 n = a->i[i+1] - diag[i] - 1; 1034 PLogFlops(2*n-1); 1035 idx = a->j + diag[i] + (!shift); 1036 v = a->a + diag[i] + (!shift); 1037 sum = xb[i]; 1038 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1039 x[i] = omega*(sum/d); 1040 } 1041 } 1042 its--; 1043 } 1044 while (its--) { 1045 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1046 for ( i=0; i<m; i++ ) { 1047 d = fshift + a->a[diag[i]+shift]; 1048 n = a->i[i+1] - a->i[i]; 1049 PLogFlops(2*n-1); 1050 idx = a->j + a->i[i] + shift; 1051 v = a->a + a->i[i] + shift; 1052 sum = b[i]; 1053 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1054 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1055 } 1056 } 1057 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1058 for ( i=m-1; i>=0; i-- ) { 1059 d = fshift + a->a[diag[i] + shift]; 1060 n = a->i[i+1] - a->i[i]; 1061 PLogFlops(2*n-1); 1062 idx = a->j + a->i[i] + shift; 1063 v = a->a + a->i[i] + shift; 1064 sum = b[i]; 1065 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1066 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1067 } 1068 } 1069 } 1070 PetscFunctionReturn(0); 1071 } 1072 1073 #undef __FUNC__ 1074 #define __FUNC__ "MatGetInfo_SeqAIJ" 1075 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1076 { 1077 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1078 1079 PetscFunctionBegin; 1080 info->rows_global = (double)a->m; 1081 info->columns_global = (double)a->n; 1082 info->rows_local = (double)a->m; 1083 info->columns_local = (double)a->n; 1084 info->block_size = 1.0; 1085 info->nz_allocated = (double)a->maxnz; 1086 info->nz_used = (double)a->nz; 1087 info->nz_unneeded = (double)(a->maxnz - a->nz); 1088 info->assemblies = (double)A->num_ass; 1089 info->mallocs = (double)a->reallocs; 1090 info->memory = A->mem; 1091 if (A->factor) { 1092 info->fill_ratio_given = A->info.fill_ratio_given; 1093 info->fill_ratio_needed = A->info.fill_ratio_needed; 1094 info->factor_mallocs = A->info.factor_mallocs; 1095 } else { 1096 info->fill_ratio_given = 0; 1097 info->fill_ratio_needed = 0; 1098 info->factor_mallocs = 0; 1099 } 1100 PetscFunctionReturn(0); 1101 } 1102 1103 extern int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,double,Mat*); 1104 extern int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1105 extern int MatLUFactor_SeqAIJ(Mat,IS,IS,double); 1106 extern int MatSolve_SeqAIJ(Mat,Vec,Vec); 1107 extern int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1108 extern int MatSolveTrans_SeqAIJ(Mat,Vec,Vec); 1109 extern int MatSolveTransAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1110 1111 #undef __FUNC__ 1112 #define __FUNC__ "MatZeroRows_SeqAIJ" 1113 int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 1114 { 1115 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1116 int i,ierr,N, *rows,m = a->m - 1,shift = a->indexshift; 1117 1118 PetscFunctionBegin; 1119 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 1120 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 1121 if (diag) { 1122 for ( i=0; i<N; i++ ) { 1123 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1124 if (a->ilen[rows[i]] > 0) { /* in case row was completely empty */ 1125 a->ilen[rows[i]] = 1; 1126 a->a[a->i[rows[i]]+shift] = *diag; 1127 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1128 } else { 1129 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1130 } 1131 } 1132 } else { 1133 for ( i=0; i<N; i++ ) { 1134 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1135 a->ilen[rows[i]] = 0; 1136 } 1137 } 1138 ISRestoreIndices(is,&rows); 1139 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1140 PetscFunctionReturn(0); 1141 } 1142 1143 #undef __FUNC__ 1144 #define __FUNC__ "MatGetSize_SeqAIJ" 1145 int MatGetSize_SeqAIJ(Mat A,int *m,int *n) 1146 { 1147 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1148 1149 PetscFunctionBegin; 1150 if (m) *m = a->m; 1151 if (n) *n = a->n; 1152 PetscFunctionReturn(0); 1153 } 1154 1155 #undef __FUNC__ 1156 #define __FUNC__ "MatGetOwnershipRange_SeqAIJ" 1157 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1158 { 1159 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1160 1161 PetscFunctionBegin; 1162 *m = 0; *n = a->m; 1163 PetscFunctionReturn(0); 1164 } 1165 1166 #undef __FUNC__ 1167 #define __FUNC__ "MatGetRow_SeqAIJ" 1168 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1169 { 1170 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1171 int *itmp,i,shift = a->indexshift; 1172 1173 PetscFunctionBegin; 1174 if (row < 0 || row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row out of range"); 1175 1176 *nz = a->i[row+1] - a->i[row]; 1177 if (v) *v = a->a + a->i[row] + shift; 1178 if (idx) { 1179 itmp = a->j + a->i[row] + shift; 1180 if (*nz && shift) { 1181 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 1182 for ( i=0; i<(*nz); i++ ) {(*idx)[i] = itmp[i] + shift;} 1183 } else if (*nz) { 1184 *idx = itmp; 1185 } 1186 else *idx = 0; 1187 } 1188 PetscFunctionReturn(0); 1189 } 1190 1191 #undef __FUNC__ 1192 #define __FUNC__ "MatRestoreRow_SeqAIJ" 1193 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1194 { 1195 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1196 1197 PetscFunctionBegin; 1198 if (idx) {if (*idx && a->indexshift) PetscFree(*idx);} 1199 PetscFunctionReturn(0); 1200 } 1201 1202 #undef __FUNC__ 1203 #define __FUNC__ "MatNorm_SeqAIJ" 1204 int MatNorm_SeqAIJ(Mat A,NormType type,double *norm) 1205 { 1206 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1207 Scalar *v = a->a; 1208 double sum = 0.0; 1209 int i, j,shift = a->indexshift; 1210 1211 PetscFunctionBegin; 1212 if (type == NORM_FROBENIUS) { 1213 for (i=0; i<a->nz; i++ ) { 1214 #if defined(USE_PETSC_COMPLEX) 1215 sum += PetscReal(PetscConj(*v)*(*v)); v++; 1216 #else 1217 sum += (*v)*(*v); v++; 1218 #endif 1219 } 1220 *norm = sqrt(sum); 1221 } else if (type == NORM_1) { 1222 double *tmp; 1223 int *jj = a->j; 1224 tmp = (double *) PetscMalloc( (a->n+1)*sizeof(double) ); CHKPTRQ(tmp); 1225 PetscMemzero(tmp,a->n*sizeof(double)); 1226 *norm = 0.0; 1227 for ( j=0; j<a->nz; j++ ) { 1228 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1229 } 1230 for ( j=0; j<a->n; j++ ) { 1231 if (tmp[j] > *norm) *norm = tmp[j]; 1232 } 1233 PetscFree(tmp); 1234 } else if (type == NORM_INFINITY) { 1235 *norm = 0.0; 1236 for ( j=0; j<a->m; j++ ) { 1237 v = a->a + a->i[j] + shift; 1238 sum = 0.0; 1239 for ( i=0; i<a->i[j+1]-a->i[j]; i++ ) { 1240 sum += PetscAbsScalar(*v); v++; 1241 } 1242 if (sum > *norm) *norm = sum; 1243 } 1244 } else { 1245 SETERRQ(PETSC_ERR_SUP,0,"No support for two norm"); 1246 } 1247 PetscFunctionReturn(0); 1248 } 1249 1250 #undef __FUNC__ 1251 #define __FUNC__ "MatTranspose_SeqAIJ" 1252 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1253 { 1254 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1255 Mat C; 1256 int i, ierr, *aj = a->j, *ai = a->i, m = a->m, len, *col; 1257 int shift = a->indexshift; 1258 Scalar *array = a->a; 1259 1260 PetscFunctionBegin; 1261 if (B == PETSC_NULL && m != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Square matrix only for in-place"); 1262 col = (int *) PetscMalloc((1+a->n)*sizeof(int)); CHKPTRQ(col); 1263 PetscMemzero(col,(1+a->n)*sizeof(int)); 1264 if (shift) { 1265 for ( i=0; i<ai[m]-1; i++ ) aj[i] -= 1; 1266 } 1267 for ( i=0; i<ai[m]+shift; i++ ) col[aj[i]] += 1; 1268 ierr = MatCreateSeqAIJ(A->comm,a->n,m,0,col,&C); CHKERRQ(ierr); 1269 PetscFree(col); 1270 for ( i=0; i<m; i++ ) { 1271 len = ai[i+1]-ai[i]; 1272 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES); CHKERRQ(ierr); 1273 array += len; aj += len; 1274 } 1275 if (shift) { 1276 for ( i=0; i<ai[m]-1; i++ ) aj[i] += 1; 1277 } 1278 1279 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1280 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1281 1282 if (B != PETSC_NULL) { 1283 *B = C; 1284 } else { 1285 PetscOps *Abops; 1286 MatOps Aops; 1287 1288 /* This isn't really an in-place transpose */ 1289 PetscFree(a->a); 1290 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 1291 if (a->diag) PetscFree(a->diag); 1292 if (a->ilen) PetscFree(a->ilen); 1293 if (a->imax) PetscFree(a->imax); 1294 if (a->solve_work) PetscFree(a->solve_work); 1295 if (a->inode.size) PetscFree(a->inode.size); 1296 PetscFree(a); 1297 1298 1299 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 1300 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 1301 1302 /* 1303 This is horrible, horrible code. We need to keep the 1304 the bops and ops Structures, copy everything from C 1305 including the function pointers.. 1306 */ 1307 PetscMemcpy(A->ops,C->ops,sizeof(struct _MatOps)); 1308 PetscMemcpy(A->bops,C->bops,sizeof(PetscOps)); 1309 Abops = A->bops; 1310 Aops = A->ops; 1311 PetscMemcpy(A,C,sizeof(struct _p_Mat)); 1312 A->bops = Abops; 1313 A->ops = Aops; 1314 A->qlist = 0; 1315 1316 PetscHeaderDestroy(C); 1317 } 1318 PetscFunctionReturn(0); 1319 } 1320 1321 #undef __FUNC__ 1322 #define __FUNC__ "MatDiagonalScale_SeqAIJ" 1323 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1324 { 1325 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1326 Scalar *l,*r,x,*v; 1327 int ierr,i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1328 1329 PetscFunctionBegin; 1330 if (ll) { 1331 /* The local size is used so that VecMPI can be passed to this routine 1332 by MatDiagonalScale_MPIAIJ */ 1333 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1334 if (m != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vector wrong length"); 1335 ierr = VecGetArray(ll,&l); CHKERRQ(ierr); 1336 v = a->a; 1337 for ( i=0; i<m; i++ ) { 1338 x = l[i]; 1339 M = a->i[i+1] - a->i[i]; 1340 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1341 } 1342 ierr = VecRestoreArray(ll,&l); CHKERRQ(ierr); 1343 PLogFlops(nz); 1344 } 1345 if (rr) { 1346 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1347 if (n != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vector wrong length"); 1348 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1349 v = a->a; jj = a->j; 1350 for ( i=0; i<nz; i++ ) { 1351 (*v++) *= r[*jj++ + shift]; 1352 } 1353 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1354 PLogFlops(nz); 1355 } 1356 PetscFunctionReturn(0); 1357 } 1358 1359 #undef __FUNC__ 1360 #define __FUNC__ "MatGetSubMatrix_SeqAIJ" 1361 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatGetSubMatrixCall scall,Mat *B) 1362 { 1363 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1364 int *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1365 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1366 register int sum,lensi; 1367 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1368 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1369 Scalar *a_new,*mat_a; 1370 Mat C; 1371 PetscTruth stride; 1372 1373 PetscFunctionBegin; 1374 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1375 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"ISrow is not sorted"); 1376 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1377 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"IScol is not sorted"); 1378 1379 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1380 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1381 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1382 1383 ierr = ISStrideGetInfo(iscol,&first,&step); CHKERRQ(ierr); 1384 ierr = ISStride(iscol,&stride); CHKERRQ(ierr); 1385 if (stride && step == 1) { 1386 /* special case of contiguous rows */ 1387 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1388 starts = lens + ncols; 1389 /* loop over new rows determining lens and starting points */ 1390 for (i=0; i<nrows; i++) { 1391 kstart = ai[irow[i]]+shift; 1392 kend = kstart + ailen[irow[i]]; 1393 for ( k=kstart; k<kend; k++ ) { 1394 if (aj[k]+shift >= first) { 1395 starts[i] = k; 1396 break; 1397 } 1398 } 1399 sum = 0; 1400 while (k < kend) { 1401 if (aj[k++]+shift >= first+ncols) break; 1402 sum++; 1403 } 1404 lens[i] = sum; 1405 } 1406 /* create submatrix */ 1407 if (scall == MAT_REUSE_MATRIX) { 1408 int n_cols,n_rows; 1409 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1410 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Reused submatrix wrong size"); 1411 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1412 C = *B; 1413 } else { 1414 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1415 } 1416 c = (Mat_SeqAIJ*) C->data; 1417 1418 /* loop over rows inserting into submatrix */ 1419 a_new = c->a; 1420 j_new = c->j; 1421 i_new = c->i; 1422 i_new[0] = -shift; 1423 for (i=0; i<nrows; i++) { 1424 ii = starts[i]; 1425 lensi = lens[i]; 1426 for ( k=0; k<lensi; k++ ) { 1427 *j_new++ = aj[ii+k] - first; 1428 } 1429 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1430 a_new += lensi; 1431 i_new[i+1] = i_new[i] + lensi; 1432 c->ilen[i] = lensi; 1433 } 1434 PetscFree(lens); 1435 } else { 1436 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1437 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1438 ssmap = smap + shift; 1439 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1440 PetscMemzero(smap,oldcols*sizeof(int)); 1441 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1442 /* determine lens of each row */ 1443 for (i=0; i<nrows; i++) { 1444 kstart = ai[irow[i]]+shift; 1445 kend = kstart + a->ilen[irow[i]]; 1446 lens[i] = 0; 1447 for ( k=kstart; k<kend; k++ ) { 1448 if (ssmap[aj[k]]) { 1449 lens[i]++; 1450 } 1451 } 1452 } 1453 /* Create and fill new matrix */ 1454 if (scall == MAT_REUSE_MATRIX) { 1455 c = (Mat_SeqAIJ *)((*B)->data); 1456 if (c->m != nrows || c->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong size"); 1457 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1458 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong no of nonzeros"); 1459 } 1460 PetscMemzero(c->ilen,c->m*sizeof(int)); 1461 C = *B; 1462 } else { 1463 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1464 } 1465 c = (Mat_SeqAIJ *)(C->data); 1466 for (i=0; i<nrows; i++) { 1467 row = irow[i]; 1468 kstart = ai[row]+shift; 1469 kend = kstart + a->ilen[row]; 1470 mat_i = c->i[i]+shift; 1471 mat_j = c->j + mat_i; 1472 mat_a = c->a + mat_i; 1473 mat_ilen = c->ilen + i; 1474 for ( k=kstart; k<kend; k++ ) { 1475 if ((tcol=ssmap[a->j[k]])) { 1476 *mat_j++ = tcol - (!shift); 1477 *mat_a++ = a->a[k]; 1478 (*mat_ilen)++; 1479 1480 } 1481 } 1482 } 1483 /* Free work space */ 1484 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1485 PetscFree(smap); PetscFree(lens); 1486 } 1487 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1488 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1489 1490 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1491 *B = C; 1492 PetscFunctionReturn(0); 1493 } 1494 1495 /* 1496 note: This can only work for identity for row and col. It would 1497 be good to check this and otherwise generate an error. 1498 */ 1499 #undef __FUNC__ 1500 #define __FUNC__ "MatILUFactor_SeqAIJ" 1501 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1502 { 1503 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1504 int ierr; 1505 Mat outA; 1506 1507 PetscFunctionBegin; 1508 if (fill != 0) SETERRQ(PETSC_ERR_SUP,0,"Only fill=0 supported"); 1509 1510 outA = inA; 1511 inA->factor = FACTOR_LU; 1512 a->row = row; 1513 a->col = col; 1514 1515 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1516 ierr = ISInvertPermutation(col,&(a->icol)); CHKERRQ(ierr); 1517 PLogObjectParent(inA,a->icol); 1518 1519 if (!a->solve_work) { /* this matrix may have been factored before */ 1520 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1521 } 1522 1523 if (!a->diag) { 1524 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1525 } 1526 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1527 PetscFunctionReturn(0); 1528 } 1529 1530 #include "pinclude/blaslapack.h" 1531 #undef __FUNC__ 1532 #define __FUNC__ "MatScale_SeqAIJ" 1533 int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1534 { 1535 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1536 int one = 1; 1537 1538 PetscFunctionBegin; 1539 BLscal_( &a->nz, alpha, a->a, &one ); 1540 PLogFlops(a->nz); 1541 PetscFunctionReturn(0); 1542 } 1543 1544 #undef __FUNC__ 1545 #define __FUNC__ "MatGetSubMatrices_SeqAIJ" 1546 int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall,Mat **B) 1547 { 1548 int ierr,i; 1549 1550 PetscFunctionBegin; 1551 if (scall == MAT_INITIAL_MATRIX) { 1552 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1553 } 1554 1555 for ( i=0; i<n; i++ ) { 1556 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1557 } 1558 PetscFunctionReturn(0); 1559 } 1560 1561 #undef __FUNC__ 1562 #define __FUNC__ "MatGetBlockSize_SeqAIJ" 1563 int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1564 { 1565 PetscFunctionBegin; 1566 *bs = 1; 1567 PetscFunctionReturn(0); 1568 } 1569 1570 #undef __FUNC__ 1571 #define __FUNC__ "MatIncreaseOverlap_SeqAIJ" 1572 int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1573 { 1574 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1575 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1576 int start, end, *ai, *aj; 1577 BT table; 1578 1579 PetscFunctionBegin; 1580 shift = a->indexshift; 1581 m = a->m; 1582 ai = a->i; 1583 aj = a->j+shift; 1584 1585 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"illegal overlap value used"); 1586 1587 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1588 ierr = BTCreate(m,table); CHKERRQ(ierr); 1589 1590 for ( i=0; i<is_max; i++ ) { 1591 /* Initialize the two local arrays */ 1592 isz = 0; 1593 BTMemzero(m,table); 1594 1595 /* Extract the indices, assume there can be duplicate entries */ 1596 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1597 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1598 1599 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1600 for ( j=0; j<n ; ++j){ 1601 if(!BTLookupSet(table, idx[j])) { nidx[isz++] = idx[j];} 1602 } 1603 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1604 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1605 1606 k = 0; 1607 for ( j=0; j<ov; j++){ /* for each overlap */ 1608 n = isz; 1609 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1610 row = nidx[k]; 1611 start = ai[row]; 1612 end = ai[row+1]; 1613 for ( l = start; l<end ; l++){ 1614 val = aj[l] + shift; 1615 if (!BTLookupSet(table,val)) {nidx[isz++] = val;} 1616 } 1617 } 1618 } 1619 ierr = ISCreateGeneral(PETSC_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1620 } 1621 BTDestroy(table); 1622 PetscFree(nidx); 1623 PetscFunctionReturn(0); 1624 } 1625 1626 /* -------------------------------------------------------------- */ 1627 #undef __FUNC__ 1628 #define __FUNC__ "MatPermute_SeqAIJ" 1629 int MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 1630 { 1631 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1632 Scalar *vwork; 1633 int i, ierr, nz, m = a->m, n = a->n, *cwork; 1634 int *row,*col,*cnew,j,*lens; 1635 IS icolp,irowp; 1636 1637 PetscFunctionBegin; 1638 ierr = ISInvertPermutation(rowp,&irowp); CHKERRQ(ierr); 1639 ierr = ISGetIndices(irowp,&row); CHKERRQ(ierr); 1640 ierr = ISInvertPermutation(colp,&icolp); CHKERRQ(ierr); 1641 ierr = ISGetIndices(icolp,&col); CHKERRQ(ierr); 1642 1643 /* determine lengths of permuted rows */ 1644 lens = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(lens); 1645 for (i=0; i<m; i++ ) { 1646 lens[row[i]] = a->i[i+1] - a->i[i]; 1647 } 1648 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1649 PetscFree(lens); 1650 1651 cnew = (int *) PetscMalloc( n*sizeof(int) ); CHKPTRQ(cnew); 1652 for (i=0; i<m; i++) { 1653 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1654 for (j=0; j<nz; j++ ) { cnew[j] = col[cwork[j]];} 1655 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES); CHKERRQ(ierr); 1656 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1657 } 1658 PetscFree(cnew); 1659 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1660 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1661 ierr = ISRestoreIndices(irowp,&row); CHKERRQ(ierr); 1662 ierr = ISRestoreIndices(icolp,&col); CHKERRQ(ierr); 1663 ierr = ISDestroy(irowp); CHKERRQ(ierr); 1664 ierr = ISDestroy(icolp); CHKERRQ(ierr); 1665 PetscFunctionReturn(0); 1666 } 1667 1668 #undef __FUNC__ 1669 #define __FUNC__ "MatPrintHelp_SeqAIJ" 1670 int MatPrintHelp_SeqAIJ(Mat A) 1671 { 1672 static int called = 0; 1673 MPI_Comm comm = A->comm; 1674 1675 PetscFunctionBegin; 1676 if (called) {PetscFunctionReturn(0);} else called = 1; 1677 (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1678 (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n"); 1679 (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n"); 1680 (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n"); 1681 (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n"); 1682 #if defined(HAVE_ESSL) 1683 (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n"); 1684 #endif 1685 PetscFunctionReturn(0); 1686 } 1687 extern int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1688 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1689 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1690 1691 /* -------------------------------------------------------------------*/ 1692 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 1693 MatGetRow_SeqAIJ,MatRestoreRow_SeqAIJ, 1694 MatMult_SeqAIJ,MatMultAdd_SeqAIJ, 1695 MatMultTrans_SeqAIJ,MatMultTransAdd_SeqAIJ, 1696 MatSolve_SeqAIJ,MatSolveAdd_SeqAIJ, 1697 MatSolveTrans_SeqAIJ,MatSolveTransAdd_SeqAIJ, 1698 MatLUFactor_SeqAIJ,0, 1699 MatRelax_SeqAIJ, 1700 MatTranspose_SeqAIJ, 1701 MatGetInfo_SeqAIJ,MatEqual_SeqAIJ, 1702 MatGetDiagonal_SeqAIJ,MatDiagonalScale_SeqAIJ,MatNorm_SeqAIJ, 1703 0,MatAssemblyEnd_SeqAIJ, 1704 MatCompress_SeqAIJ, 1705 MatSetOption_SeqAIJ,MatZeroEntries_SeqAIJ,MatZeroRows_SeqAIJ, 1706 MatLUFactorSymbolic_SeqAIJ,MatLUFactorNumeric_SeqAIJ,0,0, 1707 MatGetSize_SeqAIJ,MatGetSize_SeqAIJ,MatGetOwnershipRange_SeqAIJ, 1708 MatILUFactorSymbolic_SeqAIJ,0, 1709 0,0, 1710 MatConvertSameType_SeqAIJ,0,0, 1711 MatILUFactor_SeqAIJ,0,0, 1712 MatGetSubMatrices_SeqAIJ,MatIncreaseOverlap_SeqAIJ, 1713 MatGetValues_SeqAIJ,0, 1714 MatPrintHelp_SeqAIJ, 1715 MatScale_SeqAIJ,0,0, 1716 MatILUDTFactor_SeqAIJ, 1717 MatGetBlockSize_SeqAIJ, 1718 MatGetRowIJ_SeqAIJ, 1719 MatRestoreRowIJ_SeqAIJ, 1720 MatGetColumnIJ_SeqAIJ, 1721 MatRestoreColumnIJ_SeqAIJ, 1722 MatFDColoringCreate_SeqAIJ, 1723 MatColoringPatch_SeqAIJ, 1724 0, 1725 MatPermute_SeqAIJ, 1726 0, 1727 0, 1728 0, 1729 0, 1730 MatGetMaps_Petsc}; 1731 1732 extern int MatUseSuperLU_SeqAIJ(Mat); 1733 extern int MatUseEssl_SeqAIJ(Mat); 1734 extern int MatUseDXML_SeqAIJ(Mat); 1735 1736 #undef __FUNC__ 1737 #define __FUNC__ "MatSeqAIJSetColumnIndices_SeqAIJ" 1738 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 1739 { 1740 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1741 int i,nz,n; 1742 1743 PetscFunctionBegin; 1744 if (aij->indexshift) SETERRQ(1,1,"No support with 1 based indexing"); 1745 1746 nz = aij->maxnz; 1747 n = aij->n; 1748 for (i=0; i<nz; i++) { 1749 aij->j[i] = indices[i]; 1750 } 1751 aij->nz = nz; 1752 for ( i=0; i<n; i++ ) { 1753 aij->ilen[i] = aij->imax[i]; 1754 } 1755 1756 PetscFunctionReturn(0); 1757 } 1758 1759 #undef __FUNC__ 1760 #define __FUNC__ "MatSeqAIJSetColumnIndices" 1761 /*@ 1762 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 1763 in the matrix. 1764 1765 Input Parameters: 1766 + mat - the SeqAIJ matrix 1767 - indices - the column indices 1768 1769 Notes: 1770 This can be called if you have precomputed the nonzero structure of the 1771 matrix and want to provide it to the matrix object to improve the performance 1772 of the MatSetValues() operation. 1773 1774 You MUST have set the correct numbers of nonzeros per row in the call to 1775 MatCreateSeqAIJ(). 1776 1777 MUST be called before any calls to MatSetValues(); 1778 1779 @*/ 1780 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 1781 { 1782 int ierr,(*f)(Mat,int *); 1783 1784 PetscFunctionBegin; 1785 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1786 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void **)&f); 1787 CHKERRQ(ierr); 1788 if (f) { 1789 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1790 } else { 1791 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1792 } 1793 PetscFunctionReturn(0); 1794 } 1795 1796 #undef __FUNC__ 1797 #define __FUNC__ "MatCreateSeqAIJ" 1798 /*@C 1799 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 1800 (the default parallel PETSc format). For good matrix assembly performance 1801 the user should preallocate the matrix storage by setting the parameter nz 1802 (or the array nzz). By setting these parameters accurately, performance 1803 during matrix assembly can be increased by more than a factor of 50. 1804 1805 Collective on MPI_Comm 1806 1807 Input Parameters: 1808 + comm - MPI communicator, set to PETSC_COMM_SELF 1809 . m - number of rows 1810 . n - number of columns 1811 . nz - number of nonzeros per row (same for all rows) 1812 - nzz - array containing the number of nonzeros in the various rows 1813 (possibly different for each row) or PETSC_NULL 1814 1815 Output Parameter: 1816 . A - the matrix 1817 1818 Notes: 1819 The AIJ format (also called the Yale sparse matrix format or 1820 compressed row storage), is fully compatible with standard Fortran 77 1821 storage. That is, the stored row and column indices can begin at 1822 either one (as in Fortran) or zero. See the users' manual for details. 1823 1824 Specify the preallocated storage with either nz or nnz (not both). 1825 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1826 allocation. For large problems you MUST preallocate memory or you 1827 will get TERRIBLE performance, see the users' manual chapter on matrices. 1828 1829 By default, this format uses inodes (identical nodes) when possible, to 1830 improve numerical efficiency of matrix-vector products and solves. We 1831 search for consecutive rows with the same nonzero structure, thereby 1832 reusing matrix information to achieve increased efficiency. 1833 1834 Options Database Keys: 1835 + -mat_aij_no_inode - Do not use inodes 1836 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 1837 - -mat_aij_oneindex - Internally use indexing starting at 1 1838 rather than 0. Note that when calling MatSetValues(), 1839 the user still MUST index entries starting at 0! 1840 1841 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices() 1842 @*/ 1843 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 1844 { 1845 Mat B; 1846 Mat_SeqAIJ *b; 1847 int i, len, ierr, flg,size; 1848 1849 PetscFunctionBegin; 1850 MPI_Comm_size(comm,&size); 1851 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Comm must be of size 1"); 1852 1853 *A = 0; 1854 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,comm,MatDestroy,MatView); 1855 PLogObjectCreate(B); 1856 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 1857 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 1858 PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps)); 1859 B->ops->destroy = MatDestroy_SeqAIJ; 1860 B->ops->view = MatView_SeqAIJ; 1861 B->factor = 0; 1862 B->lupivotthreshold = 1.0; 1863 B->mapping = 0; 1864 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,&flg);CHKERRQ(ierr); 1865 b->ilu_preserve_row_sums = PETSC_FALSE; 1866 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",(int*)&b->ilu_preserve_row_sums);CHKERRQ(ierr); 1867 b->row = 0; 1868 b->col = 0; 1869 b->icol = 0; 1870 b->indexshift = 0; 1871 b->reallocs = 0; 1872 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 1873 if (flg) b->indexshift = -1; 1874 1875 b->m = m; B->m = m; B->M = m; 1876 b->n = n; B->n = n; B->N = n; 1877 1878 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 1879 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 1880 1881 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 1882 if (nnz == PETSC_NULL) { 1883 if (nz == PETSC_DEFAULT) nz = 10; 1884 else if (nz <= 0) nz = 1; 1885 for ( i=0; i<m; i++ ) b->imax[i] = nz; 1886 nz = nz*m; 1887 } else { 1888 nz = 0; 1889 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1890 } 1891 1892 /* allocate the matrix space */ 1893 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 1894 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1895 b->j = (int *) (b->a + nz); 1896 PetscMemzero(b->j,nz*sizeof(int)); 1897 b->i = b->j + nz; 1898 b->singlemalloc = 1; 1899 1900 b->i[0] = -b->indexshift; 1901 for (i=1; i<m+1; i++) { 1902 b->i[i] = b->i[i-1] + b->imax[i-1]; 1903 } 1904 1905 /* b->ilen will count nonzeros in each row so far. */ 1906 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 1907 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1908 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 1909 1910 b->nz = 0; 1911 b->maxnz = nz; 1912 b->sorted = 0; 1913 b->roworiented = 1; 1914 b->nonew = 0; 1915 b->diag = 0; 1916 b->solve_work = 0; 1917 b->spptr = 0; 1918 b->inode.node_count = 0; 1919 b->inode.size = 0; 1920 b->inode.limit = 5; 1921 b->inode.max_limit = 5; 1922 B->info.nz_unneeded = (double)b->maxnz; 1923 1924 *A = B; 1925 1926 /* SuperLU is not currently supported through PETSc */ 1927 #if defined(HAVE_SUPERLU) 1928 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 1929 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 1930 #endif 1931 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 1932 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 1933 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 1934 if (flg) { 1935 if (!b->indexshift) SETERRQ( PETSC_ERR_LIB,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 1936 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 1937 } 1938 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1939 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1940 1941 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 1942 "MatSeqAIJSetColumnIndices_SeqAIJ", 1943 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 1944 PetscFunctionReturn(0); 1945 } 1946 1947 #undef __FUNC__ 1948 #define __FUNC__ "MatConvertSameType_SeqAIJ" 1949 int MatConvertSameType_SeqAIJ(Mat A,Mat *B,int cpvalues) 1950 { 1951 Mat C; 1952 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 1953 int i,len, m = a->m,shift = a->indexshift,ierr; 1954 1955 PetscFunctionBegin; 1956 *B = 0; 1957 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,A->comm,MatDestroy,MatView); 1958 PLogObjectCreate(C); 1959 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 1960 PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps)); 1961 C->ops->destroy = MatDestroy_SeqAIJ; 1962 C->ops->view = MatView_SeqAIJ; 1963 C->factor = A->factor; 1964 c->row = 0; 1965 c->col = 0; 1966 c->icol = 0; 1967 c->indexshift = shift; 1968 C->assembled = PETSC_TRUE; 1969 1970 c->m = C->m = a->m; 1971 c->n = C->n = a->n; 1972 C->M = a->m; 1973 C->N = a->n; 1974 1975 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 1976 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 1977 for ( i=0; i<m; i++ ) { 1978 c->imax[i] = a->imax[i]; 1979 c->ilen[i] = a->ilen[i]; 1980 } 1981 1982 /* allocate the matrix space */ 1983 c->singlemalloc = 1; 1984 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 1985 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1986 c->j = (int *) (c->a + a->i[m] + shift); 1987 c->i = c->j + a->i[m] + shift; 1988 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 1989 if (m > 0) { 1990 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 1991 if (cpvalues == COPY_VALUES) { 1992 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 1993 } 1994 } 1995 1996 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1997 c->sorted = a->sorted; 1998 c->roworiented = a->roworiented; 1999 c->nonew = a->nonew; 2000 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2001 2002 if (a->diag) { 2003 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 2004 PLogObjectMemory(C,(m+1)*sizeof(int)); 2005 for ( i=0; i<m; i++ ) { 2006 c->diag[i] = a->diag[i]; 2007 } 2008 } else c->diag = 0; 2009 c->inode.limit = a->inode.limit; 2010 c->inode.max_limit = a->inode.max_limit; 2011 if (a->inode.size){ 2012 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 2013 c->inode.node_count = a->inode.node_count; 2014 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 2015 } else { 2016 c->inode.size = 0; 2017 c->inode.node_count = 0; 2018 } 2019 c->nz = a->nz; 2020 c->maxnz = a->maxnz; 2021 c->solve_work = 0; 2022 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2023 2024 *B = C; 2025 ierr = PetscObjectComposeFunction((PetscObject)C,"MatSeqAIJSetColumnIndices_C", 2026 "MatSeqAIJSetColumnIndices_SeqAIJ", 2027 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2028 PetscFunctionReturn(0); 2029 } 2030 2031 #undef __FUNC__ 2032 #define __FUNC__ "MatLoad_SeqAIJ" 2033 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 2034 { 2035 Mat_SeqAIJ *a; 2036 Mat B; 2037 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 2038 MPI_Comm comm; 2039 2040 PetscFunctionBegin; 2041 PetscObjectGetComm((PetscObject) viewer,&comm); 2042 MPI_Comm_size(comm,&size); 2043 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,0,"view must have one processor"); 2044 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 2045 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 2046 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not matrix object in file"); 2047 M = header[1]; N = header[2]; nz = header[3]; 2048 2049 if (nz < 0) { 2050 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 2051 } 2052 2053 /* read in row lengths */ 2054 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 2055 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 2056 2057 /* create our matrix */ 2058 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 2059 B = *A; 2060 a = (Mat_SeqAIJ *) B->data; 2061 shift = a->indexshift; 2062 2063 /* read in column indices and adjust for Fortran indexing*/ 2064 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT); CHKERRQ(ierr); 2065 if (shift) { 2066 for ( i=0; i<nz; i++ ) { 2067 a->j[i] += 1; 2068 } 2069 } 2070 2071 /* read in nonzero values */ 2072 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR); CHKERRQ(ierr); 2073 2074 /* set matrix "i" values */ 2075 a->i[0] = -shift; 2076 for ( i=1; i<= M; i++ ) { 2077 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2078 a->ilen[i-1] = rowlengths[i-1]; 2079 } 2080 PetscFree(rowlengths); 2081 2082 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2083 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2084 PetscFunctionReturn(0); 2085 } 2086 2087 #undef __FUNC__ 2088 #define __FUNC__ "MatEqual_SeqAIJ" 2089 int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 2090 { 2091 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 2092 2093 PetscFunctionBegin; 2094 if (B->type !=MATSEQAIJ)SETERRQ(PETSC_ERR_ARG_INCOMP,0,"Matrices must be same type"); 2095 2096 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 2097 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 2098 (a->indexshift != b->indexshift)) { 2099 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2100 } 2101 2102 /* if the a->i are the same */ 2103 if (PetscMemcmp(a->i,b->i,(a->m+1)*sizeof(int))) { 2104 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2105 } 2106 2107 /* if a->j are the same */ 2108 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 2109 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2110 } 2111 2112 /* if a->a are the same */ 2113 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 2114 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2115 } 2116 *flg = PETSC_TRUE; 2117 PetscFunctionReturn(0); 2118 2119 } 2120