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