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