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