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