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