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