1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: aij.c,v 1.285 1998/10/19 22:17:52 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 "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 } else if (format == VIEWER_FORMAT_ASCII_SYMMODU) { 385 int nzd=0, fshift=1, *sptr; 386 sptr = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(sptr); 387 for ( i=0; i<m; i++ ) { 388 sptr[i] = nzd+1; 389 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 390 if (a->j[j] >= i) { 391 #if defined(USE_PETSC_COMPLEX) 392 if (PetscImaginary(a->a[j]) != 0.0 || PetscReal(a->a[j]) != 0.0) nzd++; 393 #else 394 if (a->a[j] != 0.0) nzd++; 395 #endif 396 } 397 } 398 } 399 sptr[m] = nzd+1; 400 fprintf(fd," %d %d\n\n",m,nzd); 401 for ( i=0; i<m+1; i+=6 ) { 402 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]); 403 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]); 404 else if (i+2<m) fprintf(fd," %d %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]); 405 else if (i+1<m) fprintf(fd," %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2]); 406 else if (i<m) fprintf(fd," %d %d\n",sptr[i],sptr[i+1]); 407 else fprintf(fd," %d\n",sptr[i]); 408 } 409 fprintf(fd,"\n"); 410 PetscFree(sptr); 411 for ( i=0; i<m; i++ ) { 412 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 413 if (a->j[j] >= i) fprintf(fd," %d ",a->j[j]+fshift); 414 } 415 fprintf(fd,"\n"); 416 } 417 fprintf(fd,"\n"); 418 for ( i=0; i<m; i++ ) { 419 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 420 if (a->j[j] >= i) { 421 #if defined(USE_PETSC_COMPLEX) 422 if (PetscImaginary(a->a[j]) != 0.0 || PetscReal(a->a[j]) != 0.0) 423 fprintf(fd," %18.16e %18.16e ",PetscReal(a->a[j]),PetscImaginary(a->a[j])); 424 #else 425 if (a->a[j] != 0.0) fprintf(fd," %18.16e ",a->a[j]); 426 #endif 427 } 428 } 429 fprintf(fd,"\n"); 430 } 431 } else if (format == VIEWER_FORMAT_ASCII_DENSE) { 432 int cnt = 0,jcnt; 433 Scalar value; 434 435 for ( i=0; i<m; i++ ) { 436 jcnt = 0; 437 for ( j=0; j<a->n; j++ ) { 438 if ( jcnt++ < a->i[i+1]-a->i[i+1] && j == a->j[cnt]) { 439 value = a->a[cnt++]; 440 } else { 441 value = 0.0; 442 } 443 #if defined(USE_PETSC_COMPLEX) 444 fprintf(fd," %7.5e+%7.5e i ",PetscReal(value),PetscImaginary(value)); 445 #else 446 fprintf(fd," %7.5e ",value); 447 #endif 448 } 449 fprintf(fd,"\n"); 450 } 451 } else { 452 for ( i=0; i<m; i++ ) { 453 fprintf(fd,"row %d:",i); 454 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 455 #if defined(USE_PETSC_COMPLEX) 456 if (PetscImaginary(a->a[j]) > 0.0) { 457 fprintf(fd," %d %g + %g i",a->j[j]+shift,PetscReal(a->a[j]),PetscImaginary(a->a[j])); 458 } else if (PetscImaginary(a->a[j]) < 0.0) { 459 fprintf(fd," %d %g - %g i",a->j[j]+shift,PetscReal(a->a[j]),-PetscImaginary(a->a[j])); 460 } else { 461 fprintf(fd," %d %g ",a->j[j]+shift,PetscReal(a->a[j])); 462 } 463 #else 464 fprintf(fd," %d %g ",a->j[j]+shift,a->a[j]); 465 #endif 466 } 467 fprintf(fd,"\n"); 468 } 469 } 470 fflush(fd); 471 PetscFunctionReturn(0); 472 } 473 474 #undef __FUNC__ 475 #define __FUNC__ "MatView_SeqAIJ_Draw_Zoom" 476 int MatView_SeqAIJ_Draw_Zoom(Draw draw,void *Aa) 477 { 478 Mat A = (Mat) Aa; 479 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 480 int ierr, i,j, m = a->m, shift = a->indexshift,color; 481 int format; 482 double xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 483 Viewer viewer; 484 485 PetscFunctionBegin; 486 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*) &viewer);CHKERRQ(ierr); 487 ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr); 488 489 ierr = DrawGetCoordinates(draw,&xl,&yl,&xr,&yr); CHKERRQ(ierr); 490 /* loop over matrix elements drawing boxes */ 491 492 if (format != VIEWER_FORMAT_DRAW_CONTOUR) { 493 /* Blue for negative, Cyan for zero and Red for positive */ 494 color = DRAW_BLUE; 495 for ( i=0; i<m; i++ ) { 496 y_l = m - i - 1.0; y_r = y_l + 1.0; 497 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 498 x_l = a->j[j] + shift; x_r = x_l + 1.0; 499 #if defined(USE_PETSC_COMPLEX) 500 if (PetscReal(a->a[j]) >= 0.) continue; 501 #else 502 if (a->a[j] >= 0.) continue; 503 #endif 504 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 505 } 506 } 507 color = DRAW_CYAN; 508 for ( i=0; i<m; i++ ) { 509 y_l = m - i - 1.0; y_r = y_l + 1.0; 510 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 511 x_l = a->j[j] + shift; x_r = x_l + 1.0; 512 if (a->a[j] != 0.) continue; 513 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 514 } 515 } 516 color = DRAW_RED; 517 for ( i=0; i<m; i++ ) { 518 y_l = m - i - 1.0; y_r = y_l + 1.0; 519 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 520 x_l = a->j[j] + shift; x_r = x_l + 1.0; 521 #if defined(USE_PETSC_COMPLEX) 522 if (PetscReal(a->a[j]) <= 0.) continue; 523 #else 524 if (a->a[j] <= 0.) continue; 525 #endif 526 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 527 } 528 } 529 } else { 530 /* use contour shading to indicate magnitude of values */ 531 /* first determine max of all nonzero values */ 532 int nz = a->nz,count; 533 Draw popup; 534 double scale; 535 536 for ( i=0; i<nz; i++ ) { 537 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 538 } 539 scale = (245.0 - DRAW_BASIC_COLORS)/maxv; 540 ierr = DrawGetPopup(draw,&popup); CHKERRQ(ierr); 541 ierr = DrawScalePopup(popup,0.0,maxv); CHKERRQ(ierr); 542 count = 0; 543 for ( i=0; i<m; i++ ) { 544 y_l = m - i - 1.0; y_r = y_l + 1.0; 545 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 546 x_l = a->j[j] + shift; x_r = x_l + 1.0; 547 color = DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(a->a[count])); 548 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 549 count++; 550 } 551 } 552 } 553 PetscFunctionReturn(0); 554 } 555 556 #undef __FUNC__ 557 #define __FUNC__ "MatView_SeqAIJ_Draw" 558 int MatView_SeqAIJ_Draw(Mat A,Viewer viewer) 559 { 560 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 561 int ierr; 562 Draw draw; 563 double xr,yr,xl,yl,h,w; 564 565 PetscTruth isnull; 566 567 PetscFunctionBegin; 568 ierr = ViewerDrawGetDraw(viewer,&draw); CHKERRQ(ierr); 569 ierr = DrawIsNull(draw,&isnull); CHKERRQ(ierr); 570 if (isnull) PetscFunctionReturn(0); 571 572 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 573 xr = a->n; yr = a->m; h = yr/10.0; w = xr/10.0; 574 xr += w; yr += h; xl = -w; yl = -h; 575 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 576 ierr = DrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A); CHKERRQ(ierr); 577 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 578 PetscFunctionReturn(0); 579 } 580 581 #undef __FUNC__ 582 #define __FUNC__ "MatView_SeqAIJ" 583 int MatView_SeqAIJ(Mat A,Viewer viewer) 584 { 585 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 586 ViewerType vtype; 587 int ierr; 588 589 PetscFunctionBegin; 590 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 591 if (vtype == MATLAB_VIEWER) { 592 ierr = ViewerMatlabPutSparse_Private(viewer,a->m,a->n,a->nz,a->a,a->i,a->j); CHKERRQ(ierr); 593 } else if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER){ 594 ierr = MatView_SeqAIJ_ASCII(A,viewer); CHKERRQ(ierr); 595 } else if (vtype == BINARY_FILE_VIEWER) { 596 ierr = MatView_SeqAIJ_Binary(A,viewer); CHKERRQ(ierr); 597 } else if (vtype == DRAW_VIEWER) { 598 ierr = MatView_SeqAIJ_Draw(A,viewer); CHKERRQ(ierr); 599 } else { 600 SETERRQ(1,1,"Viewer type not supported by PETSc object"); 601 } 602 PetscFunctionReturn(0); 603 } 604 605 extern int Mat_AIJ_CheckInode(Mat); 606 #undef __FUNC__ 607 #define __FUNC__ "MatAssemblyEnd_SeqAIJ" 608 int MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 609 { 610 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 611 int fshift = 0,i,j,*ai = a->i, *aj = a->j, *imax = a->imax,ierr; 612 int m = a->m, *ip, N, *ailen = a->ilen,shift = a->indexshift,rmax = 0; 613 Scalar *aa = a->a, *ap; 614 615 PetscFunctionBegin; 616 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 617 618 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 619 for ( i=1; i<m; i++ ) { 620 /* move each row back by the amount of empty slots (fshift) before it*/ 621 fshift += imax[i-1] - ailen[i-1]; 622 rmax = PetscMax(rmax,ailen[i]); 623 if (fshift) { 624 ip = aj + ai[i] + shift; ap = aa + ai[i] + shift; 625 N = ailen[i]; 626 for ( j=0; j<N; j++ ) { 627 ip[j-fshift] = ip[j]; 628 ap[j-fshift] = ap[j]; 629 } 630 } 631 ai[i] = ai[i-1] + ailen[i-1]; 632 } 633 if (m) { 634 fshift += imax[m-1] - ailen[m-1]; 635 ai[m] = ai[m-1] + ailen[m-1]; 636 } 637 /* reset ilen and imax for each row */ 638 for ( i=0; i<m; i++ ) { 639 ailen[i] = imax[i] = ai[i+1] - ai[i]; 640 } 641 a->nz = ai[m] + shift; 642 643 /* diagonals may have moved, so kill the diagonal pointers */ 644 if (fshift && a->diag) { 645 PetscFree(a->diag); 646 PLogObjectMemory(A,-(m+1)*sizeof(int)); 647 a->diag = 0; 648 } 649 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %d X %d; storage space: %d unneeded, %d used\n", 650 m,a->n,fshift,a->nz); 651 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues is %d\n", 652 a->reallocs); 653 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Most nonzeros in any row is %d\n",rmax); 654 a->reallocs = 0; 655 A->info.nz_unneeded = (double)fshift; 656 657 /* check out for identical nodes. If found, use inode functions */ 658 ierr = Mat_AIJ_CheckInode(A); CHKERRQ(ierr); 659 PetscFunctionReturn(0); 660 } 661 662 #undef __FUNC__ 663 #define __FUNC__ "MatZeroEntries_SeqAIJ" 664 int MatZeroEntries_SeqAIJ(Mat A) 665 { 666 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 667 668 PetscFunctionBegin; 669 PetscMemzero(a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar)); 670 PetscFunctionReturn(0); 671 } 672 673 #undef __FUNC__ 674 #define __FUNC__ "MatDestroy_SeqAIJ" 675 int MatDestroy_SeqAIJ(Mat A) 676 { 677 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 678 int ierr; 679 680 PetscFunctionBegin; 681 if (--A->refct > 0) PetscFunctionReturn(0); 682 683 if (A->mapping) { 684 ierr = ISLocalToGlobalMappingDestroy(A->mapping); CHKERRQ(ierr); 685 } 686 if (A->bmapping) { 687 ierr = ISLocalToGlobalMappingDestroy(A->bmapping); CHKERRQ(ierr); 688 } 689 if (A->rmap) { 690 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 691 } 692 if (A->cmap) { 693 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 694 } 695 #if defined(USE_PETSC_LOG) 696 PLogObjectState((PetscObject)A,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 697 #endif 698 PetscFree(a->a); 699 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 700 if (a->diag) PetscFree(a->diag); 701 if (a->ilen) PetscFree(a->ilen); 702 if (a->imax) PetscFree(a->imax); 703 if (a->solve_work) PetscFree(a->solve_work); 704 if (a->inode.size) PetscFree(a->inode.size); 705 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 706 if (a->saved_values) PetscFree(a->saved_values); 707 PetscFree(a); 708 709 PLogObjectDestroy(A); 710 PetscHeaderDestroy(A); 711 PetscFunctionReturn(0); 712 } 713 714 #undef __FUNC__ 715 #define __FUNC__ "MatCompress_SeqAIJ" 716 int MatCompress_SeqAIJ(Mat A) 717 { 718 PetscFunctionBegin; 719 PetscFunctionReturn(0); 720 } 721 722 #undef __FUNC__ 723 #define __FUNC__ "MatSetOption_SeqAIJ" 724 int MatSetOption_SeqAIJ(Mat A,MatOption op) 725 { 726 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 727 728 PetscFunctionBegin; 729 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 730 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 731 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 732 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 733 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 734 else if (op == MAT_NEW_NONZERO_LOCATION_ERROR) a->nonew = -1; 735 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERROR) a->nonew = -2; 736 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 737 else if (op == MAT_ROWS_SORTED || 738 op == MAT_ROWS_UNSORTED || 739 op == MAT_SYMMETRIC || 740 op == MAT_STRUCTURALLY_SYMMETRIC || 741 op == MAT_YES_NEW_DIAGONALS || 742 op == MAT_IGNORE_OFF_PROC_ENTRIES|| 743 op == MAT_USE_HASH_TABLE) 744 PLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 745 else if (op == MAT_NO_NEW_DIAGONALS) { 746 SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS"); 747 } else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 748 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 749 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 750 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 751 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 752 else SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 753 PetscFunctionReturn(0); 754 } 755 756 #undef __FUNC__ 757 #define __FUNC__ "MatGetDiagonal_SeqAIJ" 758 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 759 { 760 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 761 int i,j, n,shift = a->indexshift,ierr; 762 Scalar *x, zero = 0.0; 763 764 PetscFunctionBegin; 765 ierr = VecSet(&zero,v);CHKERRQ(ierr); 766 ierr = VecGetArray(v,&x);;CHKERRQ(ierr); 767 ierr = VecGetLocalSize(v,&n);;CHKERRQ(ierr); 768 if (n != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming matrix and vector"); 769 for ( i=0; i<a->m; i++ ) { 770 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 771 if (a->j[j]+shift == i) { 772 x[i] = a->a[j]; 773 break; 774 } 775 } 776 } 777 ierr = VecRestoreArray(v,&x);;CHKERRQ(ierr); 778 PetscFunctionReturn(0); 779 } 780 781 /* -------------------------------------------------------*/ 782 /* Should check that shapes of vectors and matrices match */ 783 /* -------------------------------------------------------*/ 784 #undef __FUNC__ 785 #define __FUNC__ "MatMultTrans_SeqAIJ" 786 int MatMultTrans_SeqAIJ(Mat A,Vec xx,Vec yy) 787 { 788 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 789 Scalar *x, *y, *v, alpha; 790 int ierr,m = a->m, n, i, *idx, shift = a->indexshift; 791 792 PetscFunctionBegin; 793 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 794 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 795 PetscMemzero(y,a->n*sizeof(Scalar)); 796 y = y + shift; /* shift for Fortran start by 1 indexing */ 797 for ( i=0; i<m; i++ ) { 798 idx = a->j + a->i[i] + shift; 799 v = a->a + a->i[i] + shift; 800 n = a->i[i+1] - a->i[i]; 801 alpha = x[i]; 802 while (n-->0) {y[*idx++] += alpha * *v++;} 803 } 804 PLogFlops(2*a->nz - a->n); 805 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 806 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 807 PetscFunctionReturn(0); 808 } 809 810 #undef __FUNC__ 811 #define __FUNC__ "MatMultTransAdd_SeqAIJ" 812 int MatMultTransAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 813 { 814 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 815 Scalar *x, *y, *v, alpha; 816 int ierr,m = a->m, n, i, *idx,shift = a->indexshift; 817 818 PetscFunctionBegin; 819 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 820 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 821 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 822 y = y + shift; /* shift for Fortran start by 1 indexing */ 823 for ( i=0; i<m; i++ ) { 824 idx = a->j + a->i[i] + shift; 825 v = a->a + a->i[i] + shift; 826 n = a->i[i+1] - a->i[i]; 827 alpha = x[i]; 828 while (n-->0) {y[*idx++] += alpha * *v++;} 829 } 830 PLogFlops(2*a->nz); 831 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 832 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 833 PetscFunctionReturn(0); 834 } 835 836 #undef __FUNC__ 837 #define __FUNC__ "MatMult_SeqAIJ" 838 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 839 { 840 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 841 Scalar *x, *y, *v, sum; 842 int ierr,m = a->m, *idx, shift = a->indexshift,*ii; 843 #if !defined(USE_FORTRAN_KERNEL_MULTAIJ) 844 int n, i, jrow,j; 845 #endif 846 847 #if defined(HAVE_PRAGMA_DISJOINT) 848 #pragma disjoint(*x,*y,*v) 849 #endif 850 851 PetscFunctionBegin; 852 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 853 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 854 x = x + shift; /* shift for Fortran start by 1 indexing */ 855 idx = a->j; 856 v = a->a; 857 ii = a->i; 858 #if defined(USE_FORTRAN_KERNEL_MULTAIJ) 859 fortranmultaij_(&m,x,ii,idx+shift,v+shift,y); 860 #else 861 v += shift; /* shift for Fortran start by 1 indexing */ 862 idx += shift; 863 for ( i=0; i<m; i++ ) { 864 jrow = ii[i]; 865 n = ii[i+1] - jrow; 866 sum = 0.0; 867 for ( j=0; j<n; j++) { 868 sum += v[jrow]*x[idx[jrow]]; jrow++; 869 } 870 y[i] = sum; 871 } 872 #endif 873 PLogFlops(2*a->nz - m); 874 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 875 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 876 PetscFunctionReturn(0); 877 } 878 879 #undef __FUNC__ 880 #define __FUNC__ "MatMultAdd_SeqAIJ" 881 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 882 { 883 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 884 Scalar *x, *y, *z, *v, sum; 885 int ierr,m = a->m, *idx, shift = a->indexshift,*ii; 886 #if !defined(USE_FORTRAN_KERNEL_MULTADDAIJ) 887 int n,i,jrow,j; 888 #endif 889 890 PetscFunctionBegin; 891 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 892 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 893 if (zz != yy) { 894 ierr = VecGetArray(zz,&z); CHKERRQ(ierr); 895 } else { 896 z = y; 897 } 898 x = x + shift; /* shift for Fortran start by 1 indexing */ 899 idx = a->j; 900 v = a->a; 901 ii = a->i; 902 #if defined(USE_FORTRAN_KERNEL_MULTADDAIJ) 903 fortranmultaddaij_(&m,x,ii,idx+shift,v+shift,y,z); 904 #else 905 v += shift; /* shift for Fortran start by 1 indexing */ 906 idx += shift; 907 for ( i=0; i<m; i++ ) { 908 jrow = ii[i]; 909 n = ii[i+1] - jrow; 910 sum = y[i]; 911 for ( j=0; j<n; j++) { 912 sum += v[jrow]*x[idx[jrow]]; jrow++; 913 } 914 z[i] = sum; 915 } 916 #endif 917 PLogFlops(2*a->nz); 918 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 919 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 920 if (zz != yy) { 921 ierr = VecRestoreArray(zz,&z); CHKERRQ(ierr); 922 } 923 PetscFunctionReturn(0); 924 } 925 926 /* 927 Adds diagonal pointers to sparse matrix structure. 928 */ 929 930 #undef __FUNC__ 931 #define __FUNC__ "MatMarkDiag_SeqAIJ" 932 int MatMarkDiag_SeqAIJ(Mat A) 933 { 934 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 935 int i,j, *diag, m = a->m,shift = a->indexshift; 936 937 PetscFunctionBegin; 938 diag = (int *) PetscMalloc( (m+1)*sizeof(int)); CHKPTRQ(diag); 939 PLogObjectMemory(A,(m+1)*sizeof(int)); 940 for ( i=0; i<a->m; i++ ) { 941 diag[i] = a->i[i+1]; 942 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 943 if (a->j[j]+shift == i) { 944 diag[i] = j - shift; 945 break; 946 } 947 } 948 } 949 a->diag = diag; 950 PetscFunctionReturn(0); 951 } 952 953 #undef __FUNC__ 954 #define __FUNC__ "MatRelax_SeqAIJ" 955 int MatRelax_SeqAIJ(Mat A,Vec bb,double omega,MatSORType flag, 956 double fshift,int its,Vec xx) 957 { 958 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 959 Scalar *x, *b, *bs, d, *xs, sum, *v = a->a,*t,scale,*ts, *xb; 960 int ierr, *idx, *diag,n = a->n, m = a->m, i, shift = a->indexshift; 961 962 PetscFunctionBegin; 963 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 964 if (xx != bb) { 965 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 966 } else { 967 b = x; 968 } 969 970 if (!a->diag) {ierr = MatMarkDiag_SeqAIJ(A);CHKERRQ(ierr);} 971 diag = a->diag; 972 xs = x + shift; /* shifted by one for index start of a or a->j*/ 973 if (flag == SOR_APPLY_UPPER) { 974 /* apply ( U + D/omega) to the vector */ 975 bs = b + shift; 976 for ( i=0; i<m; i++ ) { 977 d = fshift + a->a[diag[i] + shift]; 978 n = a->i[i+1] - diag[i] - 1; 979 PLogFlops(2*n-1); 980 idx = a->j + diag[i] + (!shift); 981 v = a->a + diag[i] + (!shift); 982 sum = b[i]*d/omega; 983 SPARSEDENSEDOT(sum,bs,v,idx,n); 984 x[i] = sum; 985 } 986 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 987 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 988 PetscFunctionReturn(0); 989 } 990 if (flag == SOR_APPLY_LOWER) { 991 SETERRQ(PETSC_ERR_SUP,0,"SOR_APPLY_LOWER is not done"); 992 } else if (flag & SOR_EISENSTAT) { 993 /* Let A = L + U + D; where L is lower trianglar, 994 U is upper triangular, E is diagonal; This routine applies 995 996 (L + E)^{-1} A (U + E)^{-1} 997 998 to a vector efficiently using Eisenstat's trick. This is for 999 the case of SSOR preconditioner, so E is D/omega where omega 1000 is the relaxation factor. 1001 */ 1002 t = (Scalar *) PetscMalloc( m*sizeof(Scalar) ); CHKPTRQ(t); 1003 scale = (2.0/omega) - 1.0; 1004 1005 /* x = (E + U)^{-1} b */ 1006 for ( i=m-1; i>=0; i-- ) { 1007 d = fshift + a->a[diag[i] + shift]; 1008 n = a->i[i+1] - diag[i] - 1; 1009 PLogFlops(2*n-1); 1010 idx = a->j + diag[i] + (!shift); 1011 v = a->a + diag[i] + (!shift); 1012 sum = b[i]; 1013 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1014 x[i] = omega*(sum/d); 1015 } 1016 1017 /* t = b - (2*E - D)x */ 1018 v = a->a; 1019 for ( i=0; i<m; i++ ) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 1020 PLogFlops(2*m); 1021 1022 1023 /* t = (E + L)^{-1}t */ 1024 ts = t + shift; /* shifted by one for index start of a or a->j*/ 1025 diag = a->diag; 1026 for ( i=0; i<m; i++ ) { 1027 d = fshift + a->a[diag[i]+shift]; 1028 n = diag[i] - a->i[i]; 1029 PLogFlops(2*n-1); 1030 idx = a->j + a->i[i] + shift; 1031 v = a->a + a->i[i] + shift; 1032 sum = t[i]; 1033 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1034 t[i] = omega*(sum/d); 1035 } 1036 1037 /* x = x + t */ 1038 for ( i=0; i<m; i++ ) { x[i] += t[i]; } 1039 PetscFree(t); 1040 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 1041 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1042 PetscFunctionReturn(0); 1043 } 1044 if (flag & SOR_ZERO_INITIAL_GUESS) { 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 = diag[i] - 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] = omega*(sum/d); 1055 } 1056 xb = x; 1057 } else xb = b; 1058 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1059 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1060 for ( i=0; i<m; i++ ) { 1061 x[i] *= a->a[diag[i]+shift]; 1062 } 1063 PLogFlops(m); 1064 } 1065 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1066 for ( i=m-1; i>=0; i-- ) { 1067 d = fshift + a->a[diag[i] + shift]; 1068 n = a->i[i+1] - diag[i] - 1; 1069 PLogFlops(2*n-1); 1070 idx = a->j + diag[i] + (!shift); 1071 v = a->a + diag[i] + (!shift); 1072 sum = xb[i]; 1073 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1074 x[i] = omega*(sum/d); 1075 } 1076 } 1077 its--; 1078 } 1079 while (its--) { 1080 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1081 for ( i=0; i<m; i++ ) { 1082 d = fshift + a->a[diag[i]+shift]; 1083 n = a->i[i+1] - a->i[i]; 1084 PLogFlops(2*n-1); 1085 idx = a->j + a->i[i] + shift; 1086 v = a->a + a->i[i] + shift; 1087 sum = b[i]; 1088 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1089 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1090 } 1091 } 1092 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1093 for ( i=m-1; i>=0; i-- ) { 1094 d = fshift + a->a[diag[i] + shift]; 1095 n = a->i[i+1] - a->i[i]; 1096 PLogFlops(2*n-1); 1097 idx = a->j + a->i[i] + shift; 1098 v = a->a + a->i[i] + shift; 1099 sum = b[i]; 1100 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1101 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1102 } 1103 } 1104 } 1105 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 1106 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1107 PetscFunctionReturn(0); 1108 } 1109 1110 #undef __FUNC__ 1111 #define __FUNC__ "MatGetInfo_SeqAIJ" 1112 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1113 { 1114 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1115 1116 PetscFunctionBegin; 1117 info->rows_global = (double)a->m; 1118 info->columns_global = (double)a->n; 1119 info->rows_local = (double)a->m; 1120 info->columns_local = (double)a->n; 1121 info->block_size = 1.0; 1122 info->nz_allocated = (double)a->maxnz; 1123 info->nz_used = (double)a->nz; 1124 info->nz_unneeded = (double)(a->maxnz - a->nz); 1125 info->assemblies = (double)A->num_ass; 1126 info->mallocs = (double)a->reallocs; 1127 info->memory = A->mem; 1128 if (A->factor) { 1129 info->fill_ratio_given = A->info.fill_ratio_given; 1130 info->fill_ratio_needed = A->info.fill_ratio_needed; 1131 info->factor_mallocs = A->info.factor_mallocs; 1132 } else { 1133 info->fill_ratio_given = 0; 1134 info->fill_ratio_needed = 0; 1135 info->factor_mallocs = 0; 1136 } 1137 PetscFunctionReturn(0); 1138 } 1139 1140 extern int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,double,Mat*); 1141 extern int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1142 extern int MatLUFactor_SeqAIJ(Mat,IS,IS,double); 1143 extern int MatSolve_SeqAIJ(Mat,Vec,Vec); 1144 extern int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1145 extern int MatSolveTrans_SeqAIJ(Mat,Vec,Vec); 1146 extern int MatSolveTransAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1147 1148 #undef __FUNC__ 1149 #define __FUNC__ "MatZeroRows_SeqAIJ" 1150 int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 1151 { 1152 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1153 int i,ierr,N, *rows,m = a->m - 1,shift = a->indexshift; 1154 1155 PetscFunctionBegin; 1156 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 1157 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 1158 if (diag) { 1159 for ( i=0; i<N; i++ ) { 1160 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1161 if (a->ilen[rows[i]] > 0) { /* in case row was completely empty */ 1162 a->ilen[rows[i]] = 1; 1163 a->a[a->i[rows[i]]+shift] = *diag; 1164 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1165 } else { 1166 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1167 } 1168 } 1169 } else { 1170 for ( i=0; i<N; i++ ) { 1171 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1172 a->ilen[rows[i]] = 0; 1173 } 1174 } 1175 ISRestoreIndices(is,&rows); 1176 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1177 PetscFunctionReturn(0); 1178 } 1179 1180 #undef __FUNC__ 1181 #define __FUNC__ "MatGetSize_SeqAIJ" 1182 int MatGetSize_SeqAIJ(Mat A,int *m,int *n) 1183 { 1184 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1185 1186 PetscFunctionBegin; 1187 if (m) *m = a->m; 1188 if (n) *n = a->n; 1189 PetscFunctionReturn(0); 1190 } 1191 1192 #undef __FUNC__ 1193 #define __FUNC__ "MatGetOwnershipRange_SeqAIJ" 1194 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1195 { 1196 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1197 1198 PetscFunctionBegin; 1199 *m = 0; *n = a->m; 1200 PetscFunctionReturn(0); 1201 } 1202 1203 #undef __FUNC__ 1204 #define __FUNC__ "MatGetRow_SeqAIJ" 1205 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1206 { 1207 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1208 int *itmp,i,shift = a->indexshift; 1209 1210 PetscFunctionBegin; 1211 if (row < 0 || row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row out of range"); 1212 1213 *nz = a->i[row+1] - a->i[row]; 1214 if (v) *v = a->a + a->i[row] + shift; 1215 if (idx) { 1216 itmp = a->j + a->i[row] + shift; 1217 if (*nz && shift) { 1218 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 1219 for ( i=0; i<(*nz); i++ ) {(*idx)[i] = itmp[i] + shift;} 1220 } else if (*nz) { 1221 *idx = itmp; 1222 } 1223 else *idx = 0; 1224 } 1225 PetscFunctionReturn(0); 1226 } 1227 1228 #undef __FUNC__ 1229 #define __FUNC__ "MatRestoreRow_SeqAIJ" 1230 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1231 { 1232 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1233 1234 PetscFunctionBegin; 1235 if (idx) {if (*idx && a->indexshift) PetscFree(*idx);} 1236 PetscFunctionReturn(0); 1237 } 1238 1239 #undef __FUNC__ 1240 #define __FUNC__ "MatNorm_SeqAIJ" 1241 int MatNorm_SeqAIJ(Mat A,NormType type,double *norm) 1242 { 1243 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1244 Scalar *v = a->a; 1245 double sum = 0.0; 1246 int i, j,shift = a->indexshift; 1247 1248 PetscFunctionBegin; 1249 if (type == NORM_FROBENIUS) { 1250 for (i=0; i<a->nz; i++ ) { 1251 #if defined(USE_PETSC_COMPLEX) 1252 sum += PetscReal(PetscConj(*v)*(*v)); v++; 1253 #else 1254 sum += (*v)*(*v); v++; 1255 #endif 1256 } 1257 *norm = sqrt(sum); 1258 } else if (type == NORM_1) { 1259 double *tmp; 1260 int *jj = a->j; 1261 tmp = (double *) PetscMalloc( (a->n+1)*sizeof(double) ); CHKPTRQ(tmp); 1262 PetscMemzero(tmp,a->n*sizeof(double)); 1263 *norm = 0.0; 1264 for ( j=0; j<a->nz; j++ ) { 1265 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1266 } 1267 for ( j=0; j<a->n; j++ ) { 1268 if (tmp[j] > *norm) *norm = tmp[j]; 1269 } 1270 PetscFree(tmp); 1271 } else if (type == NORM_INFINITY) { 1272 *norm = 0.0; 1273 for ( j=0; j<a->m; j++ ) { 1274 v = a->a + a->i[j] + shift; 1275 sum = 0.0; 1276 for ( i=0; i<a->i[j+1]-a->i[j]; i++ ) { 1277 sum += PetscAbsScalar(*v); v++; 1278 } 1279 if (sum > *norm) *norm = sum; 1280 } 1281 } else { 1282 SETERRQ(PETSC_ERR_SUP,0,"No support for two norm"); 1283 } 1284 PetscFunctionReturn(0); 1285 } 1286 1287 #undef __FUNC__ 1288 #define __FUNC__ "MatTranspose_SeqAIJ" 1289 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1290 { 1291 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1292 Mat C; 1293 int i, ierr, *aj = a->j, *ai = a->i, m = a->m, len, *col; 1294 int shift = a->indexshift; 1295 Scalar *array = a->a; 1296 1297 PetscFunctionBegin; 1298 if (B == PETSC_NULL && m != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Square matrix only for in-place"); 1299 col = (int *) PetscMalloc((1+a->n)*sizeof(int)); CHKPTRQ(col); 1300 PetscMemzero(col,(1+a->n)*sizeof(int)); 1301 if (shift) { 1302 for ( i=0; i<ai[m]-1; i++ ) aj[i] -= 1; 1303 } 1304 for ( i=0; i<ai[m]+shift; i++ ) col[aj[i]] += 1; 1305 ierr = MatCreateSeqAIJ(A->comm,a->n,m,0,col,&C); CHKERRQ(ierr); 1306 PetscFree(col); 1307 for ( i=0; i<m; i++ ) { 1308 len = ai[i+1]-ai[i]; 1309 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES); CHKERRQ(ierr); 1310 array += len; aj += len; 1311 } 1312 if (shift) { 1313 for ( i=0; i<ai[m]-1; i++ ) aj[i] += 1; 1314 } 1315 1316 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1317 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1318 1319 if (B != PETSC_NULL) { 1320 *B = C; 1321 } else { 1322 PetscOps *Abops; 1323 MatOps Aops; 1324 1325 /* This isn't really an in-place transpose */ 1326 PetscFree(a->a); 1327 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 1328 if (a->diag) PetscFree(a->diag); 1329 if (a->ilen) PetscFree(a->ilen); 1330 if (a->imax) PetscFree(a->imax); 1331 if (a->solve_work) PetscFree(a->solve_work); 1332 if (a->inode.size) PetscFree(a->inode.size); 1333 PetscFree(a); 1334 1335 1336 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 1337 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 1338 1339 /* 1340 This is horrible, horrible code. We need to keep the 1341 the bops and ops Structures, copy everything from C 1342 including the function pointers.. 1343 */ 1344 PetscMemcpy(A->ops,C->ops,sizeof(struct _MatOps)); 1345 PetscMemcpy(A->bops,C->bops,sizeof(PetscOps)); 1346 Abops = A->bops; 1347 Aops = A->ops; 1348 PetscMemcpy(A,C,sizeof(struct _p_Mat)); 1349 A->bops = Abops; 1350 A->ops = Aops; 1351 A->qlist = 0; 1352 1353 PetscHeaderDestroy(C); 1354 } 1355 PetscFunctionReturn(0); 1356 } 1357 1358 #undef __FUNC__ 1359 #define __FUNC__ "MatDiagonalScale_SeqAIJ" 1360 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1361 { 1362 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1363 Scalar *l,*r,x,*v; 1364 int ierr,i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1365 1366 PetscFunctionBegin; 1367 if (ll) { 1368 /* The local size is used so that VecMPI can be passed to this routine 1369 by MatDiagonalScale_MPIAIJ */ 1370 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1371 if (m != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vector wrong length"); 1372 ierr = VecGetArray(ll,&l); CHKERRQ(ierr); 1373 v = a->a; 1374 for ( i=0; i<m; i++ ) { 1375 x = l[i]; 1376 M = a->i[i+1] - a->i[i]; 1377 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1378 } 1379 ierr = VecRestoreArray(ll,&l); CHKERRQ(ierr); 1380 PLogFlops(nz); 1381 } 1382 if (rr) { 1383 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1384 if (n != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vector wrong length"); 1385 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1386 v = a->a; jj = a->j; 1387 for ( i=0; i<nz; i++ ) { 1388 (*v++) *= r[*jj++ + shift]; 1389 } 1390 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1391 PLogFlops(nz); 1392 } 1393 PetscFunctionReturn(0); 1394 } 1395 1396 #undef __FUNC__ 1397 #define __FUNC__ "MatGetSubMatrix_SeqAIJ" 1398 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatGetSubMatrixCall scall,Mat *B) 1399 { 1400 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1401 int *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1402 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1403 register int sum,lensi; 1404 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1405 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1406 Scalar *a_new,*mat_a; 1407 Mat C; 1408 PetscTruth stride; 1409 1410 PetscFunctionBegin; 1411 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1412 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"ISrow is not sorted"); 1413 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1414 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"IScol is not sorted"); 1415 1416 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1417 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1418 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1419 1420 ierr = ISStrideGetInfo(iscol,&first,&step); CHKERRQ(ierr); 1421 ierr = ISStride(iscol,&stride); CHKERRQ(ierr); 1422 if (stride && step == 1) { 1423 /* special case of contiguous rows */ 1424 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1425 starts = lens + ncols; 1426 /* loop over new rows determining lens and starting points */ 1427 for (i=0; i<nrows; i++) { 1428 kstart = ai[irow[i]]+shift; 1429 kend = kstart + ailen[irow[i]]; 1430 for ( k=kstart; k<kend; k++ ) { 1431 if (aj[k]+shift >= first) { 1432 starts[i] = k; 1433 break; 1434 } 1435 } 1436 sum = 0; 1437 while (k < kend) { 1438 if (aj[k++]+shift >= first+ncols) break; 1439 sum++; 1440 } 1441 lens[i] = sum; 1442 } 1443 /* create submatrix */ 1444 if (scall == MAT_REUSE_MATRIX) { 1445 int n_cols,n_rows; 1446 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1447 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Reused submatrix wrong size"); 1448 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1449 C = *B; 1450 } else { 1451 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1452 } 1453 c = (Mat_SeqAIJ*) C->data; 1454 1455 /* loop over rows inserting into submatrix */ 1456 a_new = c->a; 1457 j_new = c->j; 1458 i_new = c->i; 1459 i_new[0] = -shift; 1460 for (i=0; i<nrows; i++) { 1461 ii = starts[i]; 1462 lensi = lens[i]; 1463 for ( k=0; k<lensi; k++ ) { 1464 *j_new++ = aj[ii+k] - first; 1465 } 1466 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1467 a_new += lensi; 1468 i_new[i+1] = i_new[i] + lensi; 1469 c->ilen[i] = lensi; 1470 } 1471 PetscFree(lens); 1472 } else { 1473 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1474 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1475 ssmap = smap + shift; 1476 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1477 PetscMemzero(smap,oldcols*sizeof(int)); 1478 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1479 /* determine lens of each row */ 1480 for (i=0; i<nrows; i++) { 1481 kstart = ai[irow[i]]+shift; 1482 kend = kstart + a->ilen[irow[i]]; 1483 lens[i] = 0; 1484 for ( k=kstart; k<kend; k++ ) { 1485 if (ssmap[aj[k]]) { 1486 lens[i]++; 1487 } 1488 } 1489 } 1490 /* Create and fill new matrix */ 1491 if (scall == MAT_REUSE_MATRIX) { 1492 c = (Mat_SeqAIJ *)((*B)->data); 1493 if (c->m != nrows || c->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong size"); 1494 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1495 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong no of nonzeros"); 1496 } 1497 PetscMemzero(c->ilen,c->m*sizeof(int)); 1498 C = *B; 1499 } else { 1500 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1501 } 1502 c = (Mat_SeqAIJ *)(C->data); 1503 for (i=0; i<nrows; i++) { 1504 row = irow[i]; 1505 kstart = ai[row]+shift; 1506 kend = kstart + a->ilen[row]; 1507 mat_i = c->i[i]+shift; 1508 mat_j = c->j + mat_i; 1509 mat_a = c->a + mat_i; 1510 mat_ilen = c->ilen + i; 1511 for ( k=kstart; k<kend; k++ ) { 1512 if ((tcol=ssmap[a->j[k]])) { 1513 *mat_j++ = tcol - (!shift); 1514 *mat_a++ = a->a[k]; 1515 (*mat_ilen)++; 1516 1517 } 1518 } 1519 } 1520 /* Free work space */ 1521 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1522 PetscFree(smap); PetscFree(lens); 1523 } 1524 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1525 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1526 1527 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1528 *B = C; 1529 PetscFunctionReturn(0); 1530 } 1531 1532 /* 1533 note: This can only work for identity for row and col. It would 1534 be good to check this and otherwise generate an error. 1535 */ 1536 #undef __FUNC__ 1537 #define __FUNC__ "MatILUFactor_SeqAIJ" 1538 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1539 { 1540 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1541 int ierr; 1542 Mat outA; 1543 1544 PetscFunctionBegin; 1545 if (fill != 0) SETERRQ(PETSC_ERR_SUP,0,"Only fill=0 supported"); 1546 1547 outA = inA; 1548 inA->factor = FACTOR_LU; 1549 a->row = row; 1550 a->col = col; 1551 1552 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1553 ierr = ISInvertPermutation(col,&(a->icol)); CHKERRQ(ierr); 1554 PLogObjectParent(inA,a->icol); 1555 1556 if (!a->solve_work) { /* this matrix may have been factored before */ 1557 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1558 } 1559 1560 if (!a->diag) { 1561 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1562 } 1563 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1564 PetscFunctionReturn(0); 1565 } 1566 1567 #include "pinclude/blaslapack.h" 1568 #undef __FUNC__ 1569 #define __FUNC__ "MatScale_SeqAIJ" 1570 int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1571 { 1572 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1573 int one = 1; 1574 1575 PetscFunctionBegin; 1576 BLscal_( &a->nz, alpha, a->a, &one ); 1577 PLogFlops(a->nz); 1578 PetscFunctionReturn(0); 1579 } 1580 1581 #undef __FUNC__ 1582 #define __FUNC__ "MatGetSubMatrices_SeqAIJ" 1583 int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall,Mat **B) 1584 { 1585 int ierr,i; 1586 1587 PetscFunctionBegin; 1588 if (scall == MAT_INITIAL_MATRIX) { 1589 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1590 } 1591 1592 for ( i=0; i<n; i++ ) { 1593 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1594 } 1595 PetscFunctionReturn(0); 1596 } 1597 1598 #undef __FUNC__ 1599 #define __FUNC__ "MatGetBlockSize_SeqAIJ" 1600 int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1601 { 1602 PetscFunctionBegin; 1603 *bs = 1; 1604 PetscFunctionReturn(0); 1605 } 1606 1607 #undef __FUNC__ 1608 #define __FUNC__ "MatIncreaseOverlap_SeqAIJ" 1609 int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1610 { 1611 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1612 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1613 int start, end, *ai, *aj; 1614 BT table; 1615 1616 PetscFunctionBegin; 1617 shift = a->indexshift; 1618 m = a->m; 1619 ai = a->i; 1620 aj = a->j+shift; 1621 1622 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"illegal overlap value used"); 1623 1624 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1625 ierr = BTCreate(m,table); CHKERRQ(ierr); 1626 1627 for ( i=0; i<is_max; i++ ) { 1628 /* Initialize the two local arrays */ 1629 isz = 0; 1630 BTMemzero(m,table); 1631 1632 /* Extract the indices, assume there can be duplicate entries */ 1633 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1634 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1635 1636 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1637 for ( j=0; j<n ; ++j){ 1638 if(!BTLookupSet(table, idx[j])) { nidx[isz++] = idx[j];} 1639 } 1640 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1641 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1642 1643 k = 0; 1644 for ( j=0; j<ov; j++){ /* for each overlap */ 1645 n = isz; 1646 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1647 row = nidx[k]; 1648 start = ai[row]; 1649 end = ai[row+1]; 1650 for ( l = start; l<end ; l++){ 1651 val = aj[l] + shift; 1652 if (!BTLookupSet(table,val)) {nidx[isz++] = val;} 1653 } 1654 } 1655 } 1656 ierr = ISCreateGeneral(PETSC_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1657 } 1658 BTDestroy(table); 1659 PetscFree(nidx); 1660 PetscFunctionReturn(0); 1661 } 1662 1663 /* -------------------------------------------------------------- */ 1664 #undef __FUNC__ 1665 #define __FUNC__ "MatPermute_SeqAIJ" 1666 int MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 1667 { 1668 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1669 Scalar *vwork; 1670 int i, ierr, nz, m = a->m, n = a->n, *cwork; 1671 int *row,*col,*cnew,j,*lens; 1672 IS icolp,irowp; 1673 1674 PetscFunctionBegin; 1675 ierr = ISInvertPermutation(rowp,&irowp); CHKERRQ(ierr); 1676 ierr = ISGetIndices(irowp,&row); CHKERRQ(ierr); 1677 ierr = ISInvertPermutation(colp,&icolp); CHKERRQ(ierr); 1678 ierr = ISGetIndices(icolp,&col); CHKERRQ(ierr); 1679 1680 /* determine lengths of permuted rows */ 1681 lens = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(lens); 1682 for (i=0; i<m; i++ ) { 1683 lens[row[i]] = a->i[i+1] - a->i[i]; 1684 } 1685 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1686 PetscFree(lens); 1687 1688 cnew = (int *) PetscMalloc( n*sizeof(int) ); CHKPTRQ(cnew); 1689 for (i=0; i<m; i++) { 1690 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1691 for (j=0; j<nz; j++ ) { cnew[j] = col[cwork[j]];} 1692 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES); CHKERRQ(ierr); 1693 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1694 } 1695 PetscFree(cnew); 1696 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1697 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1698 ierr = ISRestoreIndices(irowp,&row); CHKERRQ(ierr); 1699 ierr = ISRestoreIndices(icolp,&col); CHKERRQ(ierr); 1700 ierr = ISDestroy(irowp); CHKERRQ(ierr); 1701 ierr = ISDestroy(icolp); CHKERRQ(ierr); 1702 PetscFunctionReturn(0); 1703 } 1704 1705 #undef __FUNC__ 1706 #define __FUNC__ "MatPrintHelp_SeqAIJ" 1707 int MatPrintHelp_SeqAIJ(Mat A) 1708 { 1709 static int called = 0; 1710 MPI_Comm comm = A->comm; 1711 1712 PetscFunctionBegin; 1713 if (called) {PetscFunctionReturn(0);} else called = 1; 1714 (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1715 (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n"); 1716 (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n"); 1717 (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n"); 1718 (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n"); 1719 #if defined(HAVE_ESSL) 1720 (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n"); 1721 #endif 1722 PetscFunctionReturn(0); 1723 } 1724 extern int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1725 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1726 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1727 1728 #undef __FUNC__ 1729 #define __FUNC__ "MatCopy_SeqAIJ" 1730 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1731 { 1732 int ierr; 1733 1734 PetscFunctionBegin; 1735 if (str == SAME_NONZERO_PATTERN && B->type == MATSEQAIJ) { 1736 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1737 Mat_SeqAIJ *b = (Mat_SeqAIJ *) B->data; 1738 1739 if (a->nonew != 1) SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1740 if (b->nonew != 1) SETERRQ(1,1,"Must call MatSetOption(B,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1741 if (a->i[a->m]+a->indexshift != b->i[b->m]+a->indexshift) { 1742 SETERRQ(1,1,"Number of nonzeros in two matrices are different"); 1743 } 1744 PetscMemcpy(b->a,a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar)); 1745 } else { 1746 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1747 } 1748 PetscFunctionReturn(0); 1749 } 1750 1751 /* -------------------------------------------------------------------*/ 1752 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 1753 MatGetRow_SeqAIJ, 1754 MatRestoreRow_SeqAIJ, 1755 MatMult_SeqAIJ, 1756 MatMultAdd_SeqAIJ, 1757 MatMultTrans_SeqAIJ, 1758 MatMultTransAdd_SeqAIJ, 1759 MatSolve_SeqAIJ, 1760 MatSolveAdd_SeqAIJ, 1761 MatSolveTrans_SeqAIJ, 1762 MatSolveTransAdd_SeqAIJ, 1763 MatLUFactor_SeqAIJ, 1764 0, 1765 MatRelax_SeqAIJ, 1766 MatTranspose_SeqAIJ, 1767 MatGetInfo_SeqAIJ, 1768 MatEqual_SeqAIJ, 1769 MatGetDiagonal_SeqAIJ, 1770 MatDiagonalScale_SeqAIJ, 1771 MatNorm_SeqAIJ, 1772 0, 1773 MatAssemblyEnd_SeqAIJ, 1774 MatCompress_SeqAIJ, 1775 MatSetOption_SeqAIJ, 1776 MatZeroEntries_SeqAIJ, 1777 MatZeroRows_SeqAIJ, 1778 MatLUFactorSymbolic_SeqAIJ, 1779 MatLUFactorNumeric_SeqAIJ, 1780 0, 1781 0, 1782 MatGetSize_SeqAIJ, 1783 MatGetSize_SeqAIJ, 1784 MatGetOwnershipRange_SeqAIJ, 1785 MatILUFactorSymbolic_SeqAIJ, 1786 0, 1787 0, 1788 0, 1789 MatDuplicate_SeqAIJ, 1790 0, 1791 0, 1792 MatILUFactor_SeqAIJ, 1793 0, 1794 0, 1795 MatGetSubMatrices_SeqAIJ, 1796 MatIncreaseOverlap_SeqAIJ, 1797 MatGetValues_SeqAIJ, 1798 MatCopy_SeqAIJ, 1799 MatPrintHelp_SeqAIJ, 1800 MatScale_SeqAIJ, 1801 0, 1802 0, 1803 MatILUDTFactor_SeqAIJ, 1804 MatGetBlockSize_SeqAIJ, 1805 MatGetRowIJ_SeqAIJ, 1806 MatRestoreRowIJ_SeqAIJ, 1807 MatGetColumnIJ_SeqAIJ, 1808 MatRestoreColumnIJ_SeqAIJ, 1809 MatFDColoringCreate_SeqAIJ, 1810 MatColoringPatch_SeqAIJ, 1811 0, 1812 MatPermute_SeqAIJ, 1813 0, 1814 0, 1815 0, 1816 0, 1817 MatGetMaps_Petsc}; 1818 1819 extern int MatUseSuperLU_SeqAIJ(Mat); 1820 extern int MatUseEssl_SeqAIJ(Mat); 1821 extern int MatUseDXML_SeqAIJ(Mat); 1822 1823 EXTERN_C_BEGIN 1824 #undef __FUNC__ 1825 #define __FUNC__ "MatSeqAIJSetColumnIndices_SeqAIJ" 1826 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 1827 { 1828 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1829 int i,nz,n; 1830 1831 PetscFunctionBegin; 1832 if (aij->indexshift) SETERRQ(1,1,"No support with 1 based indexing"); 1833 1834 nz = aij->maxnz; 1835 n = aij->n; 1836 for (i=0; i<nz; i++) { 1837 aij->j[i] = indices[i]; 1838 } 1839 aij->nz = nz; 1840 for ( i=0; i<n; i++ ) { 1841 aij->ilen[i] = aij->imax[i]; 1842 } 1843 1844 PetscFunctionReturn(0); 1845 } 1846 EXTERN_C_END 1847 1848 #undef __FUNC__ 1849 #define __FUNC__ "MatSeqAIJSetColumnIndices" 1850 /*@ 1851 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 1852 in the matrix. 1853 1854 Input Parameters: 1855 + mat - the SeqAIJ matrix 1856 - indices - the column indices 1857 1858 Notes: 1859 This can be called if you have precomputed the nonzero structure of the 1860 matrix and want to provide it to the matrix object to improve the performance 1861 of the MatSetValues() operation. 1862 1863 You MUST have set the correct numbers of nonzeros per row in the call to 1864 MatCreateSeqAIJ(). 1865 1866 MUST be called before any calls to MatSetValues(); 1867 1868 @*/ 1869 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 1870 { 1871 int ierr,(*f)(Mat,int *); 1872 1873 PetscFunctionBegin; 1874 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1875 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void **)&f); 1876 CHKERRQ(ierr); 1877 if (f) { 1878 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1879 } else { 1880 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1881 } 1882 PetscFunctionReturn(0); 1883 } 1884 1885 /* ----------------------------------------------------------------------------------------*/ 1886 1887 EXTERN_C_BEGIN 1888 #undef __FUNC__ 1889 #define __FUNC__ "MatStoreValues_SeqAIJ" 1890 int MatStoreValues_SeqAIJ(Mat mat) 1891 { 1892 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1893 int nz = aij->i[aij->m]+aij->indexshift; 1894 1895 PetscFunctionBegin; 1896 if (aij->nonew != 1) { 1897 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1898 } 1899 1900 /* allocate space for values if not already there */ 1901 if (!aij->saved_values) { 1902 aij->saved_values = (Scalar *) PetscMalloc(nz*sizeof(Scalar));CHKPTRQ(aij->saved_values); 1903 } 1904 1905 /* copy values over */ 1906 PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar)); 1907 PetscFunctionReturn(0); 1908 } 1909 EXTERN_C_END 1910 1911 #undef __FUNC__ 1912 #define __FUNC__ "MatStoreValues" 1913 /*@ 1914 MatStoreValues - Stashes a copy of the matrix values; this allows, for 1915 example, reuse of the linear part of a Jacobian, while recomputing the 1916 nonlinear portion. 1917 1918 Collect on Mat 1919 1920 Input Parameters: 1921 . mat - the matrix (currently on AIJ matrices support this option) 1922 1923 Common Usage, with SNESSolve(): 1924 $ Create Jacobian matrix 1925 $ Set linear terms into matrix 1926 $ Apply boundary conditions to matrix, at this time matrix must have 1927 $ final nonzero structure (i.e. setting the nonlinear terms and applying 1928 $ boundary conditions again will not change the nonzero structure 1929 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 1930 $ ierr = MatStoreValues(mat); 1931 $ Call SNESSetJacobian() with matrix 1932 $ In your Jacobian routine 1933 $ ierr = MatRetrieveValues(mat); 1934 $ Set nonlinear terms in matrix 1935 1936 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 1937 $ // build linear portion of Jacobian 1938 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 1939 $ ierr = MatStoreValues(mat); 1940 $ loop over nonlinear iterations 1941 $ ierr = MatRetrieveValues(mat); 1942 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 1943 $ // call MatAssemblyBegin/End() on matrix 1944 $ Solve linear system with Jacobian 1945 $ endloop 1946 1947 Notes: 1948 Matrix must already be assemblied before calling this routine 1949 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 1950 calling this routine. 1951 1952 .seealso: MatRetrieveValues() 1953 1954 @*/ 1955 int MatStoreValues(Mat mat) 1956 { 1957 int ierr,(*f)(Mat); 1958 1959 PetscFunctionBegin; 1960 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1961 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1962 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1963 1964 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void **)&f); 1965 CHKERRQ(ierr); 1966 if (f) { 1967 ierr = (*f)(mat);CHKERRQ(ierr); 1968 } else { 1969 SETERRQ(1,1,"Wrong type of matrix to store values"); 1970 } 1971 PetscFunctionReturn(0); 1972 } 1973 1974 EXTERN_C_BEGIN 1975 #undef __FUNC__ 1976 #define __FUNC__ "MatRetrieveValues_SeqAIJ" 1977 int MatRetrieveValues_SeqAIJ(Mat mat) 1978 { 1979 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1980 int nz = aij->i[aij->m]+aij->indexshift; 1981 1982 PetscFunctionBegin; 1983 if (aij->nonew != 1) { 1984 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1985 } 1986 if (!aij->saved_values) { 1987 SETERRQ(1,1,"Must call MatStoreValues(A);first"); 1988 } 1989 1990 /* copy values over */ 1991 PetscMemcpy(aij->a, aij->saved_values,nz*sizeof(Scalar)); 1992 PetscFunctionReturn(0); 1993 } 1994 EXTERN_C_END 1995 1996 #undef __FUNC__ 1997 #define __FUNC__ "MatRetrieveValues" 1998 /*@ 1999 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2000 example, reuse of the linear part of a Jacobian, while recomputing the 2001 nonlinear portion. 2002 2003 Collect on Mat 2004 2005 Input Parameters: 2006 . mat - the matrix (currently on AIJ matrices support this option) 2007 2008 .seealso: MatStoreValues() 2009 2010 @*/ 2011 int MatRetrieveValues(Mat mat) 2012 { 2013 int ierr,(*f)(Mat); 2014 2015 PetscFunctionBegin; 2016 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2017 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2018 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2019 2020 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void **)&f); 2021 CHKERRQ(ierr); 2022 if (f) { 2023 ierr = (*f)(mat);CHKERRQ(ierr); 2024 } else { 2025 SETERRQ(1,1,"Wrong type of matrix to retrieve values"); 2026 } 2027 PetscFunctionReturn(0); 2028 } 2029 2030 /* --------------------------------------------------------------------------------*/ 2031 2032 #undef __FUNC__ 2033 #define __FUNC__ "MatCreateSeqAIJ" 2034 /*@C 2035 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2036 (the default parallel PETSc format). For good matrix assembly performance 2037 the user should preallocate the matrix storage by setting the parameter nz 2038 (or the array nzz). By setting these parameters accurately, performance 2039 during matrix assembly can be increased by more than a factor of 50. 2040 2041 Collective on MPI_Comm 2042 2043 Input Parameters: 2044 + comm - MPI communicator, set to PETSC_COMM_SELF 2045 . m - number of rows 2046 . n - number of columns 2047 . nz - number of nonzeros per row (same for all rows) 2048 - nzz - array containing the number of nonzeros in the various rows 2049 (possibly different for each row) or PETSC_NULL 2050 2051 Output Parameter: 2052 . A - the matrix 2053 2054 Notes: 2055 The AIJ format (also called the Yale sparse matrix format or 2056 compressed row storage), is fully compatible with standard Fortran 77 2057 storage. That is, the stored row and column indices can begin at 2058 either one (as in Fortran) or zero. See the users' manual for details. 2059 2060 Specify the preallocated storage with either nz or nnz (not both). 2061 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2062 allocation. For large problems you MUST preallocate memory or you 2063 will get TERRIBLE performance, see the users' manual chapter on matrices. 2064 2065 By default, this format uses inodes (identical nodes) when possible, to 2066 improve numerical efficiency of matrix-vector products and solves. We 2067 search for consecutive rows with the same nonzero structure, thereby 2068 reusing matrix information to achieve increased efficiency. 2069 2070 Options Database Keys: 2071 + -mat_aij_no_inode - Do not use inodes 2072 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2073 - -mat_aij_oneindex - Internally use indexing starting at 1 2074 rather than 0. Note that when calling MatSetValues(), 2075 the user still MUST index entries starting at 0! 2076 2077 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices() 2078 @*/ 2079 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 2080 { 2081 Mat B; 2082 Mat_SeqAIJ *b; 2083 int i, len, ierr, flg,size; 2084 2085 PetscFunctionBegin; 2086 MPI_Comm_size(comm,&size); 2087 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Comm must be of size 1"); 2088 2089 *A = 0; 2090 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,comm,MatDestroy,MatView); 2091 PLogObjectCreate(B); 2092 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 2093 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 2094 PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps)); 2095 B->ops->destroy = MatDestroy_SeqAIJ; 2096 B->ops->view = MatView_SeqAIJ; 2097 B->factor = 0; 2098 B->lupivotthreshold = 1.0; 2099 B->mapping = 0; 2100 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,&flg);CHKERRQ(ierr); 2101 b->ilu_preserve_row_sums = PETSC_FALSE; 2102 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",(int*)&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2103 b->row = 0; 2104 b->col = 0; 2105 b->icol = 0; 2106 b->indexshift = 0; 2107 b->reallocs = 0; 2108 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 2109 if (flg) b->indexshift = -1; 2110 2111 b->m = m; B->m = m; B->M = m; 2112 b->n = n; B->n = n; B->N = n; 2113 2114 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 2115 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 2116 2117 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 2118 if (nnz == PETSC_NULL) { 2119 if (nz == PETSC_DEFAULT) nz = 10; 2120 else if (nz <= 0) nz = 1; 2121 for ( i=0; i<m; i++ ) b->imax[i] = nz; 2122 nz = nz*m; 2123 } else { 2124 nz = 0; 2125 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 2126 } 2127 2128 /* allocate the matrix space */ 2129 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 2130 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 2131 b->j = (int *) (b->a + nz); 2132 PetscMemzero(b->j,nz*sizeof(int)); 2133 b->i = b->j + nz; 2134 b->singlemalloc = 1; 2135 2136 b->i[0] = -b->indexshift; 2137 for (i=1; i<m+1; i++) { 2138 b->i[i] = b->i[i-1] + b->imax[i-1]; 2139 } 2140 2141 /* b->ilen will count nonzeros in each row so far. */ 2142 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 2143 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2144 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 2145 2146 b->nz = 0; 2147 b->maxnz = nz; 2148 b->sorted = 0; 2149 b->roworiented = 1; 2150 b->nonew = 0; 2151 b->diag = 0; 2152 b->solve_work = 0; 2153 b->spptr = 0; 2154 b->inode.node_count = 0; 2155 b->inode.size = 0; 2156 b->inode.limit = 5; 2157 b->inode.max_limit = 5; 2158 b->saved_values = 0; 2159 B->info.nz_unneeded = (double)b->maxnz; 2160 2161 *A = B; 2162 2163 /* SuperLU is not currently supported through PETSc */ 2164 #if defined(HAVE_SUPERLU) 2165 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 2166 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 2167 #endif 2168 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 2169 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 2170 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 2171 if (flg) { 2172 if (!b->indexshift) SETERRQ( PETSC_ERR_LIB,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 2173 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 2174 } 2175 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 2176 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 2177 2178 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2179 "MatSeqAIJSetColumnIndices_SeqAIJ", 2180 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2181 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C", 2182 "MatStoreValues_SeqAIJ", 2183 (void*)MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2184 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C", 2185 "MatRetrieveValues_SeqAIJ", 2186 (void*)MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2187 PetscFunctionReturn(0); 2188 } 2189 2190 #undef __FUNC__ 2191 #define __FUNC__ "MatDuplicate_SeqAIJ" 2192 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2193 { 2194 Mat C; 2195 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 2196 int i,len, m = a->m,shift = a->indexshift,ierr; 2197 2198 PetscFunctionBegin; 2199 *B = 0; 2200 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,A->comm,MatDestroy,MatView); 2201 PLogObjectCreate(C); 2202 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 2203 PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps)); 2204 C->ops->destroy = MatDestroy_SeqAIJ; 2205 C->ops->view = MatView_SeqAIJ; 2206 C->factor = A->factor; 2207 c->row = 0; 2208 c->col = 0; 2209 c->icol = 0; 2210 c->indexshift = shift; 2211 C->assembled = PETSC_TRUE; 2212 2213 c->m = C->m = a->m; 2214 c->n = C->n = a->n; 2215 C->M = a->m; 2216 C->N = a->n; 2217 2218 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 2219 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 2220 for ( i=0; i<m; i++ ) { 2221 c->imax[i] = a->imax[i]; 2222 c->ilen[i] = a->ilen[i]; 2223 } 2224 2225 /* allocate the matrix space */ 2226 c->singlemalloc = 1; 2227 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 2228 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 2229 c->j = (int *) (c->a + a->i[m] + shift); 2230 c->i = c->j + a->i[m] + shift; 2231 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 2232 if (m > 0) { 2233 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 2234 if (cpvalues == MAT_COPY_VALUES) { 2235 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 2236 } else { 2237 PetscMemzero(c->a,(a->i[m]+shift)*sizeof(Scalar)); 2238 } 2239 } 2240 2241 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2242 c->sorted = a->sorted; 2243 c->roworiented = a->roworiented; 2244 c->nonew = a->nonew; 2245 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2246 c->saved_values = 0; 2247 2248 if (a->diag) { 2249 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 2250 PLogObjectMemory(C,(m+1)*sizeof(int)); 2251 for ( i=0; i<m; i++ ) { 2252 c->diag[i] = a->diag[i]; 2253 } 2254 } else c->diag = 0; 2255 c->inode.limit = a->inode.limit; 2256 c->inode.max_limit = a->inode.max_limit; 2257 if (a->inode.size){ 2258 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 2259 c->inode.node_count = a->inode.node_count; 2260 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 2261 } else { 2262 c->inode.size = 0; 2263 c->inode.node_count = 0; 2264 } 2265 c->nz = a->nz; 2266 c->maxnz = a->maxnz; 2267 c->solve_work = 0; 2268 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2269 2270 *B = C; 2271 ierr = PetscObjectComposeFunction((PetscObject)C,"MatSeqAIJSetColumnIndices_C", 2272 "MatSeqAIJSetColumnIndices_SeqAIJ", 2273 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2274 PetscFunctionReturn(0); 2275 } 2276 2277 #undef __FUNC__ 2278 #define __FUNC__ "MatLoad_SeqAIJ" 2279 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 2280 { 2281 Mat_SeqAIJ *a; 2282 Mat B; 2283 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 2284 MPI_Comm comm; 2285 2286 PetscFunctionBegin; 2287 PetscObjectGetComm((PetscObject) viewer,&comm); 2288 MPI_Comm_size(comm,&size); 2289 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,0,"view must have one processor"); 2290 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 2291 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 2292 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not matrix object in file"); 2293 M = header[1]; N = header[2]; nz = header[3]; 2294 2295 if (nz < 0) { 2296 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 2297 } 2298 2299 /* read in row lengths */ 2300 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 2301 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 2302 2303 /* create our matrix */ 2304 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 2305 B = *A; 2306 a = (Mat_SeqAIJ *) B->data; 2307 shift = a->indexshift; 2308 2309 /* read in column indices and adjust for Fortran indexing*/ 2310 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT); CHKERRQ(ierr); 2311 if (shift) { 2312 for ( i=0; i<nz; i++ ) { 2313 a->j[i] += 1; 2314 } 2315 } 2316 2317 /* read in nonzero values */ 2318 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR); CHKERRQ(ierr); 2319 2320 /* set matrix "i" values */ 2321 a->i[0] = -shift; 2322 for ( i=1; i<= M; i++ ) { 2323 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2324 a->ilen[i-1] = rowlengths[i-1]; 2325 } 2326 PetscFree(rowlengths); 2327 2328 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2329 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2330 PetscFunctionReturn(0); 2331 } 2332 2333 #undef __FUNC__ 2334 #define __FUNC__ "MatEqual_SeqAIJ" 2335 int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 2336 { 2337 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 2338 2339 PetscFunctionBegin; 2340 if (B->type !=MATSEQAIJ)SETERRQ(PETSC_ERR_ARG_INCOMP,0,"Matrices must be same type"); 2341 2342 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 2343 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 2344 (a->indexshift != b->indexshift)) { 2345 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2346 } 2347 2348 /* if the a->i are the same */ 2349 if (PetscMemcmp(a->i,b->i,(a->m+1)*sizeof(int))) { 2350 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2351 } 2352 2353 /* if a->j are the same */ 2354 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 2355 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2356 } 2357 2358 /* if a->a are the same */ 2359 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 2360 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2361 } 2362 *flg = PETSC_TRUE; 2363 PetscFunctionReturn(0); 2364 2365 } 2366