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