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