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