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