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