1 /*$Id: aij.c,v 1.335 1999/11/24 21:53:47 bsmith Exp bsmith $*/ 2 /* 3 Defines the basic matrix operations for the AIJ (compressed row) 4 matrix storage format. 5 */ 6 7 #include "sys.h" 8 #include "src/mat/impls/aij/seq/aij.h" 9 #include "src/vec/vecimpl.h" 10 #include "src/inline/spops.h" 11 #include "src/inline/dot.h" 12 #include "bitarray.h" 13 14 15 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 16 17 #undef __FUNC__ 18 #define __FUNC__ "MatGetRowIJ_SeqAIJ" 19 int MatGetRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *m,int **ia,int **ja,PetscTruth *done) 20 { 21 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 22 int ierr,i,ishift; 23 24 PetscFunctionBegin; 25 *m = A->m; 26 if (!ia) PetscFunctionReturn(0); 27 ishift = a->indexshift; 28 if (symmetric) { 29 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 30 } else if (oshift == 0 && ishift == -1) { 31 int nz = a->i[a->m]; 32 /* malloc space and subtract 1 from i and j indices */ 33 *ia = (int*)PetscMalloc((a->m+1)*sizeof(int));CHKPTRQ(*ia); 34 *ja = (int*)PetscMalloc((nz+1)*sizeof(int));CHKPTRQ(*ja); 35 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] - 1; 36 for (i=0; i<a->m+1; i++) (*ia)[i] = a->i[i] - 1; 37 } else if (oshift == 1 && ishift == 0) { 38 int nz = a->i[a->m] + 1; 39 /* malloc space and add 1 to i and j indices */ 40 *ia = (int*)PetscMalloc((a->m+1)*sizeof(int));CHKPTRQ(*ia); 41 *ja = (int*)PetscMalloc((nz+1)*sizeof(int));CHKPTRQ(*ja); 42 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 43 for (i=0; i<a->m+1; i++) (*ia)[i] = a->i[i] + 1; 44 } else { 45 *ia = a->i; *ja = a->j; 46 } 47 PetscFunctionReturn(0); 48 } 49 50 #undef __FUNC__ 51 #define __FUNC__ "MatRestoreRowIJ_SeqAIJ" 52 int MatRestoreRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia,int **ja,PetscTruth *done) 53 { 54 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 55 int ishift = a->indexshift,ierr; 56 57 PetscFunctionBegin; 58 if (!ia) PetscFunctionReturn(0); 59 if (symmetric || (oshift == 0 && ishift == -1) || (oshift == 1 && ishift == 0)) { 60 ierr = PetscFree(*ia);CHKERRQ(ierr); 61 ierr = PetscFree(*ja);CHKERRQ(ierr); 62 } 63 PetscFunctionReturn(0); 64 } 65 66 #undef __FUNC__ 67 #define __FUNC__ "MatGetColumnIJ_SeqAIJ" 68 int MatGetColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja,PetscTruth *done) 69 { 70 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 71 int ierr,i,ishift = a->indexshift,*collengths,*cia,*cja,n = A->n,m = A->m; 72 int nz = a->i[m]+ishift,row,*jj,mr,col; 73 74 PetscFunctionBegin; 75 *nn = A->n; 76 if (!ia) PetscFunctionReturn(0); 77 if (symmetric) { 78 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 79 } else { 80 collengths = (int*)PetscMalloc((n+1)*sizeof(int));CHKPTRQ(collengths); 81 ierr = PetscMemzero(collengths,n*sizeof(int));CHKERRQ(ierr); 82 cia = (int*)PetscMalloc((n+1)*sizeof(int));CHKPTRQ(cia); 83 cja = (int*)PetscMalloc((nz+1)*sizeof(int));CHKPTRQ(cja); 84 jj = a->j; 85 for (i=0; i<nz; i++) { 86 collengths[jj[i] + ishift]++; 87 } 88 cia[0] = oshift; 89 for (i=0; i<n; i++) { 90 cia[i+1] = cia[i] + collengths[i]; 91 } 92 ierr = PetscMemzero(collengths,n*sizeof(int));CHKERRQ(ierr); 93 jj = a->j; 94 for (row=0; row<m; row++) { 95 mr = a->i[row+1] - a->i[row]; 96 for (i=0; i<mr; i++) { 97 col = *jj++ + ishift; 98 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 99 } 100 } 101 ierr = PetscFree(collengths);CHKERRQ(ierr); 102 *ia = cia; *ja = cja; 103 } 104 PetscFunctionReturn(0); 105 } 106 107 #undef __FUNC__ 108 #define __FUNC__ "MatRestoreColumnIJ_SeqAIJ" 109 int MatRestoreColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia,int **ja,PetscTruth *done) 110 { 111 int ierr; 112 113 PetscFunctionBegin; 114 if (!ia) PetscFunctionReturn(0); 115 116 ierr = PetscFree(*ia);CHKERRQ(ierr); 117 ierr = PetscFree(*ja);CHKERRQ(ierr); 118 119 PetscFunctionReturn(0); 120 } 121 122 #define CHUNKSIZE 15 123 124 #undef __FUNC__ 125 #define __FUNC__ "MatSetValues_SeqAIJ" 126 int MatSetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 127 { 128 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 129 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted = a->sorted; 130 int *imax = a->imax,*ai = a->i,*ailen = a->ilen,roworiented = a->roworiented; 131 int *aj = a->j,nonew = a->nonew,shift = a->indexshift,ierr; 132 Scalar *ap,value,*aa = a->a; 133 PetscTruth ignorezeroentries = ((a->ignorezeroentries && is == ADD_VALUES) ? PETSC_TRUE:PETSC_FALSE); 134 135 PetscFunctionBegin; 136 for (k=0; k<m; k++) { /* loop over added rows */ 137 row = im[k]; 138 if (row < 0) continue; 139 #if defined(PETSC_USE_BOPT_g) 140 if (row >= a->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large: row %d max %d",row,a->m); 141 #endif 142 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 143 rmax = imax[row]; nrow = ailen[row]; 144 low = 0; 145 for (l=0; l<n; l++) { /* loop over added columns */ 146 if (in[l] < 0) continue; 147 #if defined(PETSC_USE_BOPT_g) 148 if (in[l] >= a->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large: col %d max %d",in[l],a->n); 149 #endif 150 col = in[l] - shift; 151 if (roworiented) { 152 value = v[l + k*n]; 153 } else { 154 value = v[k + l*m]; 155 } 156 if (value == 0.0 && ignorezeroentries) continue; 157 158 if (!sorted) low = 0; high = nrow; 159 while (high-low > 5) { 160 t = (low+high)/2; 161 if (rp[t] > col) high = t; 162 else low = t; 163 } 164 for (i=low; i<high; i++) { 165 if (rp[i] > col) break; 166 if (rp[i] == col) { 167 if (is == ADD_VALUES) ap[i] += value; 168 else ap[i] = value; 169 goto noinsert; 170 } 171 } 172 if (nonew == 1) goto noinsert; 173 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 174 if (nrow >= rmax) { 175 /* there is no extra room in row, therefore enlarge */ 176 int new_nz = ai[a->m] + CHUNKSIZE,len,*new_i,*new_j; 177 Scalar *new_a; 178 179 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 180 181 /* malloc new storage space */ 182 len = new_nz*(sizeof(int)+sizeof(Scalar))+(a->m+1)*sizeof(int); 183 new_a = (Scalar*)PetscMalloc(len);CHKPTRQ(new_a); 184 new_j = (int*)(new_a + new_nz); 185 new_i = new_j + new_nz; 186 187 /* copy over old data into new slots */ 188 for (ii=0; ii<row+1; ii++) {new_i[ii] = ai[ii];} 189 for (ii=row+1; ii<a->m+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 190 ierr = PetscMemcpy(new_j,aj,(ai[row]+nrow+shift)*sizeof(int));CHKERRQ(ierr); 191 len = (new_nz - CHUNKSIZE - ai[row] - nrow - shift); 192 ierr = PetscMemcpy(new_j+ai[row]+shift+nrow+CHUNKSIZE,aj+ai[row]+shift+nrow,len*sizeof(int));CHKERRQ(ierr); 193 ierr = PetscMemcpy(new_a,aa,(ai[row]+nrow+shift)*sizeof(Scalar));CHKERRQ(ierr); 194 ierr = PetscMemcpy(new_a+ai[row]+shift+nrow+CHUNKSIZE,aa+ai[row]+shift+nrow,len*sizeof(Scalar));CHKERRQ(ierr); 195 /* free up old matrix storage */ 196 ierr = PetscFree(a->a);CHKERRQ(ierr); 197 if (!a->singlemalloc) { 198 ierr = PetscFree(a->i);CHKERRQ(ierr); 199 ierr = PetscFree(a->j);CHKERRQ(ierr); 200 } 201 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 202 a->singlemalloc = PETSC_TRUE; 203 204 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 205 rmax = imax[row] = imax[row] + CHUNKSIZE; 206 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + sizeof(Scalar))); 207 a->maxnz += CHUNKSIZE; 208 a->reallocs++; 209 } 210 N = nrow++ - 1; a->nz++; 211 /* shift up all the later entries in this row */ 212 for (ii=N; ii>=i; ii--) { 213 rp[ii+1] = rp[ii]; 214 ap[ii+1] = ap[ii]; 215 } 216 rp[i] = col; 217 ap[i] = value; 218 noinsert:; 219 low = i + 1; 220 } 221 ailen[row] = nrow; 222 } 223 PetscFunctionReturn(0); 224 } 225 226 #undef __FUNC__ 227 #define __FUNC__ "MatGetValues_SeqAIJ" 228 int MatGetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 229 { 230 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 231 int *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 232 int *ai = a->i,*ailen = a->ilen,shift = a->indexshift; 233 Scalar *ap,*aa = a->a,zero = 0.0; 234 235 PetscFunctionBegin; 236 for (k=0; k<m; k++) { /* loop over rows */ 237 row = im[k]; 238 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 239 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 240 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 241 nrow = ailen[row]; 242 for (l=0; l<n; l++) { /* loop over columns */ 243 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 244 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 245 col = in[l] - shift; 246 high = nrow; low = 0; /* assume unsorted */ 247 while (high-low > 5) { 248 t = (low+high)/2; 249 if (rp[t] > col) high = t; 250 else low = t; 251 } 252 for (i=low; i<high; i++) { 253 if (rp[i] > col) break; 254 if (rp[i] == col) { 255 *v++ = ap[i]; 256 goto finished; 257 } 258 } 259 *v++ = zero; 260 finished:; 261 } 262 } 263 PetscFunctionReturn(0); 264 } 265 266 267 #undef __FUNC__ 268 #define __FUNC__ "MatView_SeqAIJ_Binary" 269 int MatView_SeqAIJ_Binary(Mat A,Viewer viewer) 270 { 271 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 272 int i,fd,*col_lens,ierr; 273 274 PetscFunctionBegin; 275 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 276 col_lens = (int*)PetscMalloc((4+a->m)*sizeof(int));CHKPTRQ(col_lens); 277 col_lens[0] = MAT_COOKIE; 278 col_lens[1] = a->m; 279 col_lens[2] = a->n; 280 col_lens[3] = a->nz; 281 282 /* store lengths of each row and write (including header) to file */ 283 for (i=0; i<a->m; i++) { 284 col_lens[4+i] = a->i[i+1] - a->i[i]; 285 } 286 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,PETSC_INT,1);CHKERRQ(ierr); 287 ierr = PetscFree(col_lens);CHKERRQ(ierr); 288 289 /* store column indices (zero start index) */ 290 if (a->indexshift) { 291 for (i=0; i<a->nz; i++) a->j[i]--; 292 } 293 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,0);CHKERRQ(ierr); 294 if (a->indexshift) { 295 for (i=0; i<a->nz; i++) a->j[i]++; 296 } 297 298 /* store nonzero values */ 299 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,0);CHKERRQ(ierr); 300 PetscFunctionReturn(0); 301 } 302 303 #undef __FUNC__ 304 #define __FUNC__ "MatView_SeqAIJ_ASCII" 305 int MatView_SeqAIJ_ASCII(Mat A,Viewer viewer) 306 { 307 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 308 int ierr,i,j,m = a->m,shift = a->indexshift,format; 309 char *outputname; 310 311 PetscFunctionBegin; 312 ierr = ViewerGetOutputname(viewer,&outputname);CHKERRQ(ierr); 313 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 314 if (format == VIEWER_FORMAT_ASCII_INFO_LONG || format == VIEWER_FORMAT_ASCII_INFO) { 315 if (a->inode.size) { 316 ierr = ViewerASCIIPrintf(viewer,"using I-node routines: found %d nodes, limit used is %d\n",a->inode.node_count,a->inode.limit);CHKERRQ(ierr); 317 } else { 318 ierr = ViewerASCIIPrintf(viewer,"not using I-node routines\n");CHKERRQ(ierr); 319 } 320 } else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 321 int nofinalvalue = 0; 322 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != a->n-!shift)) { 323 nofinalvalue = 1; 324 } 325 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 326 ierr = ViewerASCIIPrintf(viewer,"%% Size = %d %d \n",m,a->n);CHKERRQ(ierr); 327 ierr = ViewerASCIIPrintf(viewer,"%% Nonzeros = %d \n",a->nz);CHKERRQ(ierr); 328 ierr = ViewerASCIIPrintf(viewer,"zzz = zeros(%d,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 329 ierr = ViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 330 331 for (i=0; i<m; i++) { 332 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 333 #if defined(PETSC_USE_COMPLEX) 334 ierr = ViewerASCIIPrintf(viewer,"%d %d %18.16e + %18.16e i \n",i+1,a->j[j]+!shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 335 #else 336 ierr = ViewerASCIIPrintf(viewer,"%d %d %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 337 #endif 338 } 339 } 340 if (nofinalvalue) { 341 ierr = ViewerASCIIPrintf(viewer,"%d %d %18.16e\n",m,a->n,0.0);CHKERRQ(ierr); 342 } 343 if (outputname) {ierr = ViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",outputname);CHKERRQ(ierr);} 344 else {ierr = ViewerASCIIPrintf(viewer,"];\n M = spconvert(zzz);\n");CHKERRQ(ierr);} 345 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 346 } else if (format == VIEWER_FORMAT_ASCII_COMMON) { 347 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 348 for (i=0; i<m; i++) { 349 ierr = ViewerASCIIPrintf(viewer,"row %d:",i);CHKERRQ(ierr); 350 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 351 #if defined(PETSC_USE_COMPLEX) 352 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 353 ierr = ViewerASCIIPrintf(viewer," %d %g + %g i",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 354 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 355 ierr = ViewerASCIIPrintf(viewer," %d %g - %g i",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 356 } else if (PetscRealPart(a->a[j]) != 0.0) { 357 ierr = ViewerASCIIPrintf(viewer," %d %g ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 358 } 359 #else 360 if (a->a[j] != 0.0) {ierr = ViewerASCIIPrintf(viewer," %d %g ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 361 #endif 362 } 363 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 364 } 365 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 366 } else if (format == VIEWER_FORMAT_ASCII_SYMMODU) { 367 int nzd=0,fshift=1,*sptr; 368 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 369 sptr = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(sptr); 370 for (i=0; i<m; i++) { 371 sptr[i] = nzd+1; 372 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 373 if (a->j[j] >= i) { 374 #if defined(PETSC_USE_COMPLEX) 375 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 376 #else 377 if (a->a[j] != 0.0) nzd++; 378 #endif 379 } 380 } 381 } 382 sptr[m] = nzd+1; 383 ierr = ViewerASCIIPrintf(viewer," %d %d\n\n",m,nzd);CHKERRQ(ierr); 384 for (i=0; i<m+1; i+=6) { 385 if (i+4<m) {ierr = ViewerASCIIPrintf(viewer," %d %d %d %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5]);CHKERRQ(ierr);} 386 else if (i+3<m) {ierr = ViewerASCIIPrintf(viewer," %d %d %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4]);CHKERRQ(ierr);} 387 else if (i+2<m) {ierr = ViewerASCIIPrintf(viewer," %d %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]);CHKERRQ(ierr);} 388 else if (i+1<m) {ierr = ViewerASCIIPrintf(viewer," %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 389 else if (i<m) {ierr = ViewerASCIIPrintf(viewer," %d %d\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 390 else {ierr = ViewerASCIIPrintf(viewer," %d\n",sptr[i]);CHKERRQ(ierr);} 391 } 392 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 393 ierr = PetscFree(sptr);CHKERRQ(ierr); 394 for (i=0; i<m; i++) { 395 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 396 if (a->j[j] >= i) {ierr = ViewerASCIIPrintf(viewer," %d ",a->j[j]+fshift);CHKERRQ(ierr);} 397 } 398 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 399 } 400 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 401 for (i=0; i<m; i++) { 402 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 403 if (a->j[j] >= i) { 404 #if defined(PETSC_USE_COMPLEX) 405 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 406 ierr = ViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 407 } 408 #else 409 if (a->a[j] != 0.0) {ierr = ViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 410 #endif 411 } 412 } 413 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 414 } 415 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 416 } else if (format == VIEWER_FORMAT_ASCII_DENSE) { 417 int cnt = 0,jcnt; 418 Scalar value; 419 420 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 421 for (i=0; i<m; i++) { 422 jcnt = 0; 423 for (j=0; j<a->n; j++) { 424 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 425 value = a->a[cnt++]; 426 jcnt++; 427 } else { 428 value = 0.0; 429 } 430 #if defined(PETSC_USE_COMPLEX) 431 ierr = ViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 432 #else 433 ierr = ViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 434 #endif 435 } 436 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 437 } 438 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 439 } else { 440 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 441 for (i=0; i<m; i++) { 442 ierr = ViewerASCIIPrintf(viewer,"row %d:",i);CHKERRQ(ierr); 443 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 444 #if defined(PETSC_USE_COMPLEX) 445 if (PetscImaginaryPart(a->a[j]) > 0.0) { 446 ierr = ViewerASCIIPrintf(viewer," %d %g + %g i",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 447 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 448 ierr = ViewerASCIIPrintf(viewer," %d %g - %g i",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 449 } else { 450 ierr = ViewerASCIIPrintf(viewer," %d %g ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 451 } 452 #else 453 ierr = ViewerASCIIPrintf(viewer," %d %g ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 454 #endif 455 } 456 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 457 } 458 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 459 } 460 ierr = ViewerFlush(viewer);CHKERRQ(ierr); 461 PetscFunctionReturn(0); 462 } 463 464 #undef __FUNC__ 465 #define __FUNC__ "MatView_SeqAIJ_Draw_Zoom" 466 int MatView_SeqAIJ_Draw_Zoom(Draw draw,void *Aa) 467 { 468 Mat A = (Mat) Aa; 469 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 470 int ierr,i,j,m = a->m,shift = a->indexshift,color,rank; 471 int format; 472 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 473 Viewer viewer; 474 MPI_Comm comm; 475 476 PetscFunctionBegin; 477 /* 478 This is nasty. If this is called from an originally parallel matrix 479 then all processes call this,but only the first has the matrix so the 480 rest should return immediately. 481 */ 482 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 483 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 484 if (rank) PetscFunctionReturn(0); 485 486 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 487 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 488 489 ierr = DrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 490 /* loop over matrix elements drawing boxes */ 491 492 if (format != VIEWER_FORMAT_DRAW_CONTOUR) { 493 /* Blue for negative, Cyan for zero and Red for positive */ 494 color = DRAW_BLUE; 495 for (i=0; i<m; i++) { 496 y_l = m - i - 1.0; y_r = y_l + 1.0; 497 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 498 x_l = a->j[j] + shift; x_r = x_l + 1.0; 499 #if defined(PETSC_USE_COMPLEX) 500 if (PetscRealPart(a->a[j]) >= 0.) continue; 501 #else 502 if (a->a[j] >= 0.) continue; 503 #endif 504 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 505 } 506 } 507 color = DRAW_CYAN; 508 for (i=0; i<m; i++) { 509 y_l = m - i - 1.0; y_r = y_l + 1.0; 510 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 511 x_l = a->j[j] + shift; x_r = x_l + 1.0; 512 if (a->a[j] != 0.) continue; 513 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 514 } 515 } 516 color = DRAW_RED; 517 for (i=0; i<m; i++) { 518 y_l = m - i - 1.0; y_r = y_l + 1.0; 519 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 520 x_l = a->j[j] + shift; x_r = x_l + 1.0; 521 #if defined(PETSC_USE_COMPLEX) 522 if (PetscRealPart(a->a[j]) <= 0.) continue; 523 #else 524 if (a->a[j] <= 0.) continue; 525 #endif 526 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 527 } 528 } 529 } else { 530 /* use contour shading to indicate magnitude of values */ 531 /* first determine max of all nonzero values */ 532 int nz = a->nz,count; 533 Draw popup; 534 PetscReal scale; 535 536 for (i=0; i<nz; i++) { 537 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 538 } 539 scale = (245.0 - DRAW_BASIC_COLORS)/maxv; 540 ierr = DrawGetPopup(draw,&popup);CHKERRQ(ierr); 541 if (popup) {ierr = DrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 542 count = 0; 543 for (i=0; i<m; i++) { 544 y_l = m - i - 1.0; y_r = y_l + 1.0; 545 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 546 x_l = a->j[j] + shift; x_r = x_l + 1.0; 547 color = DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(a->a[count])); 548 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 549 count++; 550 } 551 } 552 } 553 PetscFunctionReturn(0); 554 } 555 556 #undef __FUNC__ 557 #define __FUNC__ "MatView_SeqAIJ_Draw" 558 int MatView_SeqAIJ_Draw(Mat A,Viewer viewer) 559 { 560 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 561 int ierr; 562 Draw draw; 563 PetscReal xr,yr,xl,yl,h,w; 564 PetscTruth isnull; 565 566 PetscFunctionBegin; 567 ierr = ViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 568 ierr = DrawIsNull(draw,&isnull);CHKERRQ(ierr); 569 if (isnull) PetscFunctionReturn(0); 570 571 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 572 xr = a->n; yr = a->m; h = yr/10.0; w = xr/10.0; 573 xr += w; yr += h; xl = -w; yl = -h; 574 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 575 ierr = DrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 576 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 577 PetscFunctionReturn(0); 578 } 579 580 #undef __FUNC__ 581 #define __FUNC__ "MatView_SeqAIJ" 582 int MatView_SeqAIJ(Mat A,Viewer viewer) 583 { 584 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 585 int ierr; 586 PetscTruth issocket,isascii,isbinary,isdraw; 587 588 PetscFunctionBegin; 589 ierr = PetscTypeCompare((PetscObject)viewer,SOCKET_VIEWER,&issocket);CHKERRQ(ierr); 590 ierr = PetscTypeCompare((PetscObject)viewer,ASCII_VIEWER,&isascii);CHKERRQ(ierr); 591 ierr = PetscTypeCompare((PetscObject)viewer,BINARY_VIEWER,&isbinary);CHKERRQ(ierr); 592 ierr = PetscTypeCompare((PetscObject)viewer,DRAW_VIEWER,&isdraw);CHKERRQ(ierr); 593 if (issocket) { 594 ierr = ViewerSocketPutSparse_Private(viewer,a->m,a->n,a->nz,a->a,a->i,a->j);CHKERRQ(ierr); 595 } else if (isascii) { 596 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 597 } else if (isbinary) { 598 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 599 } else if (isdraw) { 600 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 601 } else { 602 SETERRQ1(1,1,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 603 } 604 PetscFunctionReturn(0); 605 } 606 607 extern int Mat_AIJ_CheckInode(Mat); 608 #undef __FUNC__ 609 #define __FUNC__ "MatAssemblyEnd_SeqAIJ" 610 int MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 611 { 612 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 613 int fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax,ierr; 614 int m = a->m,*ip,N,*ailen = a->ilen,shift = a->indexshift,rmax = 0; 615 Scalar *aa = a->a,*ap; 616 617 PetscFunctionBegin; 618 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 619 620 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 621 for (i=1; i<m; i++) { 622 /* move each row back by the amount of empty slots (fshift) before it*/ 623 fshift += imax[i-1] - ailen[i-1]; 624 rmax = PetscMax(rmax,ailen[i]); 625 if (fshift) { 626 ip = aj + ai[i] + shift; 627 ap = aa + ai[i] + shift; 628 N = ailen[i]; 629 for (j=0; j<N; j++) { 630 ip[j-fshift] = ip[j]; 631 ap[j-fshift] = ap[j]; 632 } 633 } 634 ai[i] = ai[i-1] + ailen[i-1]; 635 } 636 if (m) { 637 fshift += imax[m-1] - ailen[m-1]; 638 ai[m] = ai[m-1] + ailen[m-1]; 639 } 640 /* reset ilen and imax for each row */ 641 for (i=0; i<m; i++) { 642 ailen[i] = imax[i] = ai[i+1] - ai[i]; 643 } 644 a->nz = ai[m] + shift; 645 646 /* diagonals may have moved, so kill the diagonal pointers */ 647 if (fshift && a->diag) { 648 ierr = PetscFree(a->diag);CHKERRQ(ierr); 649 PLogObjectMemory(A,-(m+1)*sizeof(int)); 650 a->diag = 0; 651 } 652 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %d X %d; storage space: %d unneeded,%d used\n",m,a->n,fshift,a->nz); 653 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues() is %d\n",a->reallocs); 654 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Most nonzeros in any row is %d\n",rmax); 655 a->reallocs = 0; 656 A->info.nz_unneeded = (double)fshift; 657 a->rmax = rmax; 658 659 /* check out for identical nodes. If found, use inode functions */ 660 ierr = Mat_AIJ_CheckInode(A);CHKERRQ(ierr); 661 PetscFunctionReturn(0); 662 } 663 664 #undef __FUNC__ 665 #define __FUNC__ "MatZeroEntries_SeqAIJ" 666 int MatZeroEntries_SeqAIJ(Mat A) 667 { 668 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 669 int ierr; 670 671 PetscFunctionBegin; 672 ierr = PetscMemzero(a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar));CHKERRQ(ierr); 673 PetscFunctionReturn(0); 674 } 675 676 #undef __FUNC__ 677 #define __FUNC__ "MatDestroy_SeqAIJ" 678 int MatDestroy_SeqAIJ(Mat A) 679 { 680 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 681 int ierr; 682 683 PetscFunctionBegin; 684 685 if (A->mapping) { 686 ierr = ISLocalToGlobalMappingDestroy(A->mapping);CHKERRQ(ierr); 687 } 688 if (A->bmapping) { 689 ierr = ISLocalToGlobalMappingDestroy(A->bmapping);CHKERRQ(ierr); 690 } 691 if (A->rmap) { 692 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 693 } 694 if (A->cmap) { 695 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 696 } 697 if (a->idiag) {ierr = PetscFree(a->idiag);CHKERRQ(ierr);} 698 #if defined(PETSC_USE_LOG) 699 PLogObjectState((PetscObject)A,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 700 #endif 701 if (a->freedata) { 702 ierr = PetscFree(a->a);CHKERRQ(ierr); 703 if (!a->singlemalloc) { 704 ierr = PetscFree(a->i);CHKERRQ(ierr); 705 ierr = PetscFree(a->j);CHKERRQ(ierr); 706 } 707 } 708 if (a->row) { 709 ierr = ISDestroy(a->row);CHKERRQ(ierr); 710 } 711 if (a->col) { 712 ierr = ISDestroy(a->col);CHKERRQ(ierr); 713 } 714 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 715 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 716 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 717 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 718 if (a->inode.size) {ierr = PetscFree(a->inode.size);CHKERRQ(ierr);} 719 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 720 if (a->saved_values) {ierr = PetscFree(a->saved_values);CHKERRQ(ierr);} 721 ierr = PetscFree(a);CHKERRQ(ierr); 722 723 PLogObjectDestroy(A); 724 PetscHeaderDestroy(A); 725 PetscFunctionReturn(0); 726 } 727 728 #undef __FUNC__ 729 #define __FUNC__ "MatCompress_SeqAIJ" 730 int MatCompress_SeqAIJ(Mat A) 731 { 732 PetscFunctionBegin; 733 PetscFunctionReturn(0); 734 } 735 736 #undef __FUNC__ 737 #define __FUNC__ "MatSetOption_SeqAIJ" 738 int MatSetOption_SeqAIJ(Mat A,MatOption op) 739 { 740 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 741 742 PetscFunctionBegin; 743 if (op == MAT_ROW_ORIENTED) a->roworiented = PETSC_TRUE; 744 else if (op == MAT_KEEP_ZEROED_ROWS) a->keepzeroedrows = PETSC_TRUE; 745 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = PETSC_FALSE; 746 else if (op == MAT_COLUMNS_SORTED) a->sorted = PETSC_TRUE; 747 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = PETSC_FALSE; 748 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 749 else if (op == MAT_NEW_NONZERO_LOCATION_ERR) a->nonew = -1; 750 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERR) a->nonew = -2; 751 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 752 else if (op == MAT_IGNORE_ZERO_ENTRIES) a->ignorezeroentries = PETSC_TRUE; 753 else if (op == MAT_ROWS_SORTED || 754 op == MAT_ROWS_UNSORTED || 755 op == MAT_SYMMETRIC || 756 op == MAT_STRUCTURALLY_SYMMETRIC || 757 op == MAT_YES_NEW_DIAGONALS || 758 op == MAT_IGNORE_OFF_PROC_ENTRIES|| 759 op == MAT_USE_HASH_TABLE) 760 PLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 761 else if (op == MAT_NO_NEW_DIAGONALS) { 762 SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS"); 763 } else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 764 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 765 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 766 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 767 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 768 else SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 769 PetscFunctionReturn(0); 770 } 771 772 #undef __FUNC__ 773 #define __FUNC__ "MatGetDiagonal_SeqAIJ" 774 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 775 { 776 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 777 int i,j,n,shift = a->indexshift,ierr; 778 Scalar *x,zero = 0.0; 779 780 PetscFunctionBegin; 781 ierr = VecSet(&zero,v);CHKERRQ(ierr); 782 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 783 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 784 if (n != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming matrix and vector"); 785 for (i=0; i<a->m; i++) { 786 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 787 if (a->j[j]+shift == i) { 788 x[i] = a->a[j]; 789 break; 790 } 791 } 792 } 793 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 794 PetscFunctionReturn(0); 795 } 796 797 /* -------------------------------------------------------*/ 798 /* Should check that shapes of vectors and matrices match */ 799 /* -------------------------------------------------------*/ 800 #undef __FUNC__ 801 #define __FUNC__ "MatMultTranspose_SeqAIJ" 802 int MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 803 { 804 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 805 Scalar *x,*y,*v,alpha,zero = 0.0; 806 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 807 808 PetscFunctionBegin; 809 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 810 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 811 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 812 y = y + shift; /* shift for Fortran start by 1 indexing */ 813 for (i=0; i<m; i++) { 814 idx = a->j + a->i[i] + shift; 815 v = a->a + a->i[i] + shift; 816 n = a->i[i+1] - a->i[i]; 817 alpha = x[i]; 818 while (n-->0) {y[*idx++] += alpha * *v++;} 819 } 820 PLogFlops(2*a->nz - a->n); 821 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 822 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 823 PetscFunctionReturn(0); 824 } 825 826 #undef __FUNC__ 827 #define __FUNC__ "MatMultTransposeAdd_SeqAIJ" 828 int MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 829 { 830 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 831 Scalar *x,*y,*v,alpha; 832 int ierr,m = a->m,n,i,*idx,shift = a->indexshift; 833 834 PetscFunctionBegin; 835 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 836 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 837 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 838 y = y + shift; /* shift for Fortran start by 1 indexing */ 839 for (i=0; i<m; i++) { 840 idx = a->j + a->i[i] + shift; 841 v = a->a + a->i[i] + shift; 842 n = a->i[i+1] - a->i[i]; 843 alpha = x[i]; 844 while (n-->0) {y[*idx++] += alpha * *v++;} 845 } 846 PLogFlops(2*a->nz); 847 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 848 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 849 PetscFunctionReturn(0); 850 } 851 852 #undef __FUNC__ 853 #define __FUNC__ "MatMult_SeqAIJ" 854 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 855 { 856 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 857 Scalar *x,*y,*v,sum; 858 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 859 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 860 int n,i,jrow,j; 861 #endif 862 863 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 864 #pragma disjoint(*x,*y,*v) 865 #endif 866 867 PetscFunctionBegin; 868 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 869 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 870 x = x + shift; /* shift for Fortran start by 1 indexing */ 871 idx = a->j; 872 v = a->a; 873 ii = a->i; 874 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 875 fortranmultaij_(&m,x,ii,idx+shift,v+shift,y); 876 #else 877 v += shift; /* shift for Fortran start by 1 indexing */ 878 idx += shift; 879 for (i=0; i<m; i++) { 880 jrow = ii[i]; 881 n = ii[i+1] - jrow; 882 sum = 0.0; 883 for (j=0; j<n; j++) { 884 sum += v[jrow]*x[idx[jrow]]; jrow++; 885 } 886 y[i] = sum; 887 } 888 #endif 889 PLogFlops(2*a->nz - m); 890 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 891 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 892 PetscFunctionReturn(0); 893 } 894 895 #undef __FUNC__ 896 #define __FUNC__ "MatMultAdd_SeqAIJ" 897 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 898 { 899 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 900 Scalar *x,*y,*z,*v,sum; 901 int ierr,m = a->m,*idx,shift = a->indexshift,*ii; 902 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 903 int n,i,jrow,j; 904 #endif 905 906 PetscFunctionBegin; 907 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 908 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 909 if (zz != yy) { 910 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 911 } else { 912 z = y; 913 } 914 x = x + shift; /* shift for Fortran start by 1 indexing */ 915 idx = a->j; 916 v = a->a; 917 ii = a->i; 918 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 919 fortranmultaddaij_(&m,x,ii,idx+shift,v+shift,y,z); 920 #else 921 v += shift; /* shift for Fortran start by 1 indexing */ 922 idx += shift; 923 for (i=0; i<m; i++) { 924 jrow = ii[i]; 925 n = ii[i+1] - jrow; 926 sum = y[i]; 927 for (j=0; j<n; j++) { 928 sum += v[jrow]*x[idx[jrow]]; jrow++; 929 } 930 z[i] = sum; 931 } 932 #endif 933 PLogFlops(2*a->nz); 934 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 935 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 936 if (zz != yy) { 937 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 938 } 939 PetscFunctionReturn(0); 940 } 941 942 /* 943 Adds diagonal pointers to sparse matrix structure. 944 */ 945 #undef __FUNC__ 946 #define __FUNC__ "MatMarkDiagonal_SeqAIJ" 947 int MatMarkDiagonal_SeqAIJ(Mat A) 948 { 949 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 950 int i,j,*diag,m = a->m,shift = a->indexshift; 951 952 PetscFunctionBegin; 953 if (a->diag) PetscFunctionReturn(0); 954 955 diag = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(diag); 956 PLogObjectMemory(A,(m+1)*sizeof(int)); 957 for (i=0; i<a->m; i++) { 958 diag[i] = a->i[i+1]; 959 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 960 if (a->j[j]+shift == i) { 961 diag[i] = j - shift; 962 break; 963 } 964 } 965 } 966 a->diag = diag; 967 PetscFunctionReturn(0); 968 } 969 970 /* 971 Checks for missing diagonals 972 */ 973 #undef __FUNC__ 974 #define __FUNC__ "MatMissingDiagonal_SeqAIJ" 975 int MatMissingDiagonal_SeqAIJ(Mat A) 976 { 977 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 978 int *diag,*jj = a->j,i,shift = a->indexshift,ierr; 979 980 PetscFunctionBegin; 981 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 982 diag = a->diag; 983 for (i=0; i<a->m; i++) { 984 if (jj[diag[i]+shift] != i-shift) { 985 SETERRQ1(1,1,"Matrix is missing diagonal number %d",i); 986 } 987 } 988 PetscFunctionReturn(0); 989 } 990 991 #undef __FUNC__ 992 #define __FUNC__ "MatRelax_SeqAIJ" 993 int MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,int its,Vec xx) 994 { 995 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 996 Scalar *x,*b,*bs, d,*xs,sum,*v = a->a,*t=0,scale,*ts,*xb,*idiag=0; 997 int ierr,*idx,*diag,n = a->n,m = a->m,i,shift = a->indexshift; 998 999 PetscFunctionBegin; 1000 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1001 if (xx != bb) { 1002 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 1003 } else { 1004 b = x; 1005 } 1006 1007 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 1008 diag = a->diag; 1009 xs = x + shift; /* shifted by one for index start of a or a->j*/ 1010 if (flag == SOR_APPLY_UPPER) { 1011 /* apply (U + D/omega) to the vector */ 1012 bs = b + shift; 1013 for (i=0; i<m; i++) { 1014 d = fshift + a->a[diag[i] + shift]; 1015 n = a->i[i+1] - diag[i] - 1; 1016 PLogFlops(2*n-1); 1017 idx = a->j + diag[i] + (!shift); 1018 v = a->a + diag[i] + (!shift); 1019 sum = b[i]*d/omega; 1020 SPARSEDENSEDOT(sum,bs,v,idx,n); 1021 x[i] = sum; 1022 } 1023 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1024 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1025 PetscFunctionReturn(0); 1026 } 1027 1028 /* setup workspace for Eisenstat */ 1029 if (flag & SOR_EISENSTAT) { 1030 if (!a->idiag) { 1031 a->idiag = (Scalar*)PetscMalloc(2*m*sizeof(Scalar));CHKPTRQ(a->idiag); 1032 a->ssor = a->idiag + m; 1033 v = a->a; 1034 for (i=0; i<m; i++) { a->idiag[i] = 1.0/v[diag[i]];} 1035 } 1036 t = a->ssor; 1037 idiag = a->idiag; 1038 } 1039 /* Let A = L + U + D; where L is lower trianglar, 1040 U is upper triangular, E is diagonal; This routine applies 1041 1042 (L + E)^{-1} A (U + E)^{-1} 1043 1044 to a vector efficiently using Eisenstat's trick. This is for 1045 the case of SSOR preconditioner, so E is D/omega where omega 1046 is the relaxation factor. 1047 */ 1048 1049 if (flag == SOR_APPLY_LOWER) { 1050 SETERRQ(PETSC_ERR_SUP,0,"SOR_APPLY_LOWER is not done"); 1051 } else if ((flag & SOR_EISENSTAT) && omega == 1.0 && shift == 0 && fshift == 0.0) { 1052 /* special case for omega = 1.0 saves flops and some integer ops */ 1053 Scalar *v2; 1054 1055 v2 = a->a; 1056 /* x = (E + U)^{-1} b */ 1057 for (i=m-1; i>=0; i--) { 1058 n = a->i[i+1] - diag[i] - 1; 1059 idx = a->j + diag[i] + 1; 1060 v = a->a + diag[i] + 1; 1061 sum = b[i]; 1062 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1063 x[i] = sum*idiag[i]; 1064 1065 /* t = b - (2*E - D)x */ 1066 t[i] = b[i] - (v2[diag[i]])*x[i]; 1067 } 1068 1069 /* t = (E + L)^{-1}t */ 1070 diag = a->diag; 1071 for (i=0; i<m; i++) { 1072 n = diag[i] - a->i[i]; 1073 idx = a->j + a->i[i]; 1074 v = a->a + a->i[i]; 1075 sum = t[i]; 1076 SPARSEDENSEMDOT(sum,t,v,idx,n); 1077 t[i] = sum*idiag[i]; 1078 1079 /* x = x + t */ 1080 x[i] += t[i]; 1081 } 1082 1083 PLogFlops(3*m-1 + 2*a->nz); 1084 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1085 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1086 PetscFunctionReturn(0); 1087 } else if (flag & SOR_EISENSTAT) { 1088 /* Let A = L + U + D; where L is lower trianglar, 1089 U is upper triangular, E is diagonal; This routine applies 1090 1091 (L + E)^{-1} A (U + E)^{-1} 1092 1093 to a vector efficiently using Eisenstat's trick. This is for 1094 the case of SSOR preconditioner, so E is D/omega where omega 1095 is the relaxation factor. 1096 */ 1097 scale = (2.0/omega) - 1.0; 1098 1099 /* x = (E + U)^{-1} b */ 1100 for (i=m-1; i>=0; i--) { 1101 d = fshift + a->a[diag[i] + shift]; 1102 n = a->i[i+1] - diag[i] - 1; 1103 idx = a->j + diag[i] + (!shift); 1104 v = a->a + diag[i] + (!shift); 1105 sum = b[i]; 1106 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1107 x[i] = omega*(sum/d); 1108 } 1109 1110 /* t = b - (2*E - D)x */ 1111 v = a->a; 1112 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 1113 1114 /* t = (E + L)^{-1}t */ 1115 ts = t + shift; /* shifted by one for index start of a or a->j*/ 1116 diag = a->diag; 1117 for (i=0; i<m; i++) { 1118 d = fshift + a->a[diag[i]+shift]; 1119 n = diag[i] - a->i[i]; 1120 idx = a->j + a->i[i] + shift; 1121 v = a->a + a->i[i] + shift; 1122 sum = t[i]; 1123 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1124 t[i] = omega*(sum/d); 1125 /* x = x + t */ 1126 x[i] += t[i]; 1127 } 1128 1129 PLogFlops(6*m-1 + 2*a->nz); 1130 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1131 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1132 PetscFunctionReturn(0); 1133 } 1134 if (flag & SOR_ZERO_INITIAL_GUESS) { 1135 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1136 for (i=0; i<m; i++) { 1137 d = fshift + a->a[diag[i]+shift]; 1138 n = diag[i] - a->i[i]; 1139 PLogFlops(2*n-1); 1140 idx = a->j + a->i[i] + shift; 1141 v = a->a + a->i[i] + shift; 1142 sum = b[i]; 1143 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1144 x[i] = omega*(sum/d); 1145 } 1146 xb = x; 1147 } else xb = b; 1148 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1149 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1150 for (i=0; i<m; i++) { 1151 x[i] *= a->a[diag[i]+shift]; 1152 } 1153 PLogFlops(m); 1154 } 1155 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1156 for (i=m-1; i>=0; i--) { 1157 d = fshift + a->a[diag[i] + shift]; 1158 n = a->i[i+1] - diag[i] - 1; 1159 PLogFlops(2*n-1); 1160 idx = a->j + diag[i] + (!shift); 1161 v = a->a + diag[i] + (!shift); 1162 sum = xb[i]; 1163 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1164 x[i] = omega*(sum/d); 1165 } 1166 } 1167 its--; 1168 } 1169 while (its--) { 1170 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1171 for (i=0; i<m; i++) { 1172 d = fshift + a->a[diag[i]+shift]; 1173 n = a->i[i+1] - a->i[i]; 1174 PLogFlops(2*n-1); 1175 idx = a->j + a->i[i] + shift; 1176 v = a->a + a->i[i] + shift; 1177 sum = b[i]; 1178 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1179 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1180 } 1181 } 1182 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1183 for (i=m-1; i>=0; i--) { 1184 d = fshift + a->a[diag[i] + shift]; 1185 n = a->i[i+1] - a->i[i]; 1186 PLogFlops(2*n-1); 1187 idx = a->j + a->i[i] + shift; 1188 v = a->a + a->i[i] + shift; 1189 sum = b[i]; 1190 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1191 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1192 } 1193 } 1194 } 1195 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1196 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1197 PetscFunctionReturn(0); 1198 } 1199 1200 #undef __FUNC__ 1201 #define __FUNC__ "MatGetInfo_SeqAIJ" 1202 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1203 { 1204 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1205 1206 PetscFunctionBegin; 1207 info->rows_global = (double)a->m; 1208 info->columns_global = (double)a->n; 1209 info->rows_local = (double)a->m; 1210 info->columns_local = (double)a->n; 1211 info->block_size = 1.0; 1212 info->nz_allocated = (double)a->maxnz; 1213 info->nz_used = (double)a->nz; 1214 info->nz_unneeded = (double)(a->maxnz - a->nz); 1215 info->assemblies = (double)A->num_ass; 1216 info->mallocs = (double)a->reallocs; 1217 info->memory = A->mem; 1218 if (A->factor) { 1219 info->fill_ratio_given = A->info.fill_ratio_given; 1220 info->fill_ratio_needed = A->info.fill_ratio_needed; 1221 info->factor_mallocs = A->info.factor_mallocs; 1222 } else { 1223 info->fill_ratio_given = 0; 1224 info->fill_ratio_needed = 0; 1225 info->factor_mallocs = 0; 1226 } 1227 PetscFunctionReturn(0); 1228 } 1229 1230 extern int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,PetscReal,Mat*); 1231 extern int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1232 extern int MatLUFactor_SeqAIJ(Mat,IS,IS,PetscReal); 1233 extern int MatSolve_SeqAIJ(Mat,Vec,Vec); 1234 extern int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1235 extern int MatSolveTranspose_SeqAIJ(Mat,Vec,Vec); 1236 extern int MatSolveTransposeAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1237 1238 #undef __FUNC__ 1239 #define __FUNC__ "MatZeroRows_SeqAIJ" 1240 int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 1241 { 1242 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1243 int i,ierr,N,*rows,m = a->m - 1,shift = a->indexshift; 1244 1245 PetscFunctionBegin; 1246 ierr = ISGetSize(is,&N);CHKERRQ(ierr); 1247 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1248 if (a->keepzeroedrows) { 1249 for (i=0; i<N; i++) { 1250 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1251 ierr = PetscMemzero(&a->a[a->i[rows[i]]+shift],a->ilen[rows[i]]*sizeof(Scalar));CHKERRQ(ierr); 1252 } 1253 if (diag) { 1254 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1255 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1256 for (i=0; i<N; i++) { 1257 a->a[a->diag[rows[i]]] = *diag; 1258 } 1259 } 1260 } else { 1261 if (diag) { 1262 for (i=0; i<N; i++) { 1263 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1264 if (a->ilen[rows[i]] > 0) { 1265 a->ilen[rows[i]] = 1; 1266 a->a[a->i[rows[i]]+shift] = *diag; 1267 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1268 } else { /* in case row was completely empty */ 1269 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1270 } 1271 } 1272 } else { 1273 for (i=0; i<N; i++) { 1274 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1275 a->ilen[rows[i]] = 0; 1276 } 1277 } 1278 } 1279 ierr = ISRestoreIndices(is,&rows);CHKERRQ(ierr); 1280 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1281 PetscFunctionReturn(0); 1282 } 1283 1284 #undef __FUNC__ 1285 #define __FUNC__ "MatGetSize_SeqAIJ" 1286 int MatGetSize_SeqAIJ(Mat A,int *m,int *n) 1287 { 1288 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1289 1290 PetscFunctionBegin; 1291 if (m) *m = a->m; 1292 if (n) *n = a->n; 1293 PetscFunctionReturn(0); 1294 } 1295 1296 #undef __FUNC__ 1297 #define __FUNC__ "MatGetOwnershipRange_SeqAIJ" 1298 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1299 { 1300 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1301 1302 PetscFunctionBegin; 1303 *m = 0; *n = a->m; 1304 PetscFunctionReturn(0); 1305 } 1306 1307 #undef __FUNC__ 1308 #define __FUNC__ "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) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row out of range"); 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__ "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__ "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__ "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__ "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__ "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__ "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,&(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__ "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__ "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__ "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__ "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__ "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,&irowp);CHKERRQ(ierr); 1799 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1800 ierr = ISInvertPermutation(colp,&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__ "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 PetscFunctionReturn(0); 1847 } 1848 extern int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg); 1849 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1850 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1851 extern int MatILUDTFactor_SeqAIJ(Mat,MatILUInfo*,IS,IS,Mat*); 1852 #undef __FUNC__ 1853 #define __FUNC__ "MatCopy_SeqAIJ" 1854 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1855 { 1856 int ierr; 1857 1858 PetscFunctionBegin; 1859 if (str == SAME_NONZERO_PATTERN && B->type == MATSEQAIJ) { 1860 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1861 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1862 1863 if (a->i[a->m]+a->indexshift != b->i[b->m]+a->indexshift) { 1864 SETERRQ(1,1,"Number of nonzeros in two matrices are different"); 1865 } 1866 ierr = PetscMemcpy(b->a,a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar));CHKERRQ(ierr); 1867 } else { 1868 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1869 } 1870 PetscFunctionReturn(0); 1871 } 1872 1873 /* -------------------------------------------------------------------*/ 1874 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 1875 MatGetRow_SeqAIJ, 1876 MatRestoreRow_SeqAIJ, 1877 MatMult_SeqAIJ, 1878 MatMultAdd_SeqAIJ, 1879 MatMultTranspose_SeqAIJ, 1880 MatMultTransposeAdd_SeqAIJ, 1881 MatSolve_SeqAIJ, 1882 MatSolveAdd_SeqAIJ, 1883 MatSolveTranspose_SeqAIJ, 1884 MatSolveTransposeAdd_SeqAIJ, 1885 MatLUFactor_SeqAIJ, 1886 0, 1887 MatRelax_SeqAIJ, 1888 MatTranspose_SeqAIJ, 1889 MatGetInfo_SeqAIJ, 1890 MatEqual_SeqAIJ, 1891 MatGetDiagonal_SeqAIJ, 1892 MatDiagonalScale_SeqAIJ, 1893 MatNorm_SeqAIJ, 1894 0, 1895 MatAssemblyEnd_SeqAIJ, 1896 MatCompress_SeqAIJ, 1897 MatSetOption_SeqAIJ, 1898 MatZeroEntries_SeqAIJ, 1899 MatZeroRows_SeqAIJ, 1900 MatLUFactorSymbolic_SeqAIJ, 1901 MatLUFactorNumeric_SeqAIJ, 1902 0, 1903 0, 1904 MatGetSize_SeqAIJ, 1905 MatGetSize_SeqAIJ, 1906 MatGetOwnershipRange_SeqAIJ, 1907 MatILUFactorSymbolic_SeqAIJ, 1908 0, 1909 0, 1910 0, 1911 MatDuplicate_SeqAIJ, 1912 0, 1913 0, 1914 MatILUFactor_SeqAIJ, 1915 0, 1916 0, 1917 MatGetSubMatrices_SeqAIJ, 1918 MatIncreaseOverlap_SeqAIJ, 1919 MatGetValues_SeqAIJ, 1920 MatCopy_SeqAIJ, 1921 MatPrintHelp_SeqAIJ, 1922 MatScale_SeqAIJ, 1923 0, 1924 0, 1925 MatILUDTFactor_SeqAIJ, 1926 MatGetBlockSize_SeqAIJ, 1927 MatGetRowIJ_SeqAIJ, 1928 MatRestoreRowIJ_SeqAIJ, 1929 MatGetColumnIJ_SeqAIJ, 1930 MatRestoreColumnIJ_SeqAIJ, 1931 MatFDColoringCreate_SeqAIJ, 1932 MatColoringPatch_SeqAIJ, 1933 0, 1934 MatPermute_SeqAIJ, 1935 0, 1936 0, 1937 0, 1938 0, 1939 MatGetMaps_Petsc}; 1940 1941 extern int MatUseSuperLU_SeqAIJ(Mat); 1942 extern int MatUseEssl_SeqAIJ(Mat); 1943 extern int MatUseDXML_SeqAIJ(Mat); 1944 1945 EXTERN_C_BEGIN 1946 #undef __FUNC__ 1947 #define __FUNC__ "MatSeqAIJSetColumnIndices_SeqAIJ" 1948 1949 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 1950 { 1951 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1952 int i,nz,n; 1953 1954 PetscFunctionBegin; 1955 if (aij->indexshift) SETERRQ(1,1,"No support with 1 based indexing"); 1956 1957 nz = aij->maxnz; 1958 n = aij->n; 1959 for (i=0; i<nz; i++) { 1960 aij->j[i] = indices[i]; 1961 } 1962 aij->nz = nz; 1963 for (i=0; i<n; i++) { 1964 aij->ilen[i] = aij->imax[i]; 1965 } 1966 1967 PetscFunctionReturn(0); 1968 } 1969 EXTERN_C_END 1970 1971 #undef __FUNC__ 1972 #define __FUNC__ "MatSeqAIJSetColumnIndices" 1973 /*@ 1974 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 1975 in the matrix. 1976 1977 Input Parameters: 1978 + mat - the SeqAIJ matrix 1979 - indices - the column indices 1980 1981 Level: advanced 1982 1983 Notes: 1984 This can be called if you have precomputed the nonzero structure of the 1985 matrix and want to provide it to the matrix object to improve the performance 1986 of the MatSetValues() operation. 1987 1988 You MUST have set the correct numbers of nonzeros per row in the call to 1989 MatCreateSeqAIJ(). 1990 1991 MUST be called before any calls to MatSetValues(); 1992 1993 @*/ 1994 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 1995 { 1996 int ierr,(*f)(Mat,int *); 1997 1998 PetscFunctionBegin; 1999 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2000 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void **)&f);CHKERRQ(ierr); 2001 if (f) { 2002 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2003 } else { 2004 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 2005 } 2006 PetscFunctionReturn(0); 2007 } 2008 2009 /* ----------------------------------------------------------------------------------------*/ 2010 2011 EXTERN_C_BEGIN 2012 #undef __FUNC__ 2013 #define __FUNC__ "MatStoreValues_SeqAIJ" 2014 int MatStoreValues_SeqAIJ(Mat mat) 2015 { 2016 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2017 int nz = aij->i[aij->m]+aij->indexshift,ierr; 2018 2019 PetscFunctionBegin; 2020 if (aij->nonew != 1) { 2021 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2022 } 2023 2024 /* allocate space for values if not already there */ 2025 if (!aij->saved_values) { 2026 aij->saved_values = (Scalar*)PetscMalloc(nz*sizeof(Scalar));CHKPTRQ(aij->saved_values); 2027 } 2028 2029 /* copy values over */ 2030 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar));CHKERRQ(ierr); 2031 PetscFunctionReturn(0); 2032 } 2033 EXTERN_C_END 2034 2035 #undef __FUNC__ 2036 #define __FUNC__ "MatStoreValues" 2037 /*@ 2038 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2039 example, reuse of the linear part of a Jacobian, while recomputing the 2040 nonlinear portion. 2041 2042 Collect on Mat 2043 2044 Input Parameters: 2045 . mat - the matrix (currently on AIJ matrices support this option) 2046 2047 Level: advanced 2048 2049 Common Usage, with SNESSolve(): 2050 $ Create Jacobian matrix 2051 $ Set linear terms into matrix 2052 $ Apply boundary conditions to matrix, at this time matrix must have 2053 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2054 $ boundary conditions again will not change the nonzero structure 2055 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2056 $ ierr = MatStoreValues(mat); 2057 $ Call SNESSetJacobian() with matrix 2058 $ In your Jacobian routine 2059 $ ierr = MatRetrieveValues(mat); 2060 $ Set nonlinear terms in matrix 2061 2062 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2063 $ // build linear portion of Jacobian 2064 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2065 $ ierr = MatStoreValues(mat); 2066 $ loop over nonlinear iterations 2067 $ ierr = MatRetrieveValues(mat); 2068 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2069 $ // call MatAssemblyBegin/End() on matrix 2070 $ Solve linear system with Jacobian 2071 $ endloop 2072 2073 Notes: 2074 Matrix must already be assemblied before calling this routine 2075 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2076 calling this routine. 2077 2078 .seealso: MatRetrieveValues() 2079 2080 @*/ 2081 int MatStoreValues(Mat mat) 2082 { 2083 int ierr,(*f)(Mat); 2084 2085 PetscFunctionBegin; 2086 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2087 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2088 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2089 2090 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void **)&f);CHKERRQ(ierr); 2091 if (f) { 2092 ierr = (*f)(mat);CHKERRQ(ierr); 2093 } else { 2094 SETERRQ(1,1,"Wrong type of matrix to store values"); 2095 } 2096 PetscFunctionReturn(0); 2097 } 2098 2099 EXTERN_C_BEGIN 2100 #undef __FUNC__ 2101 #define __FUNC__ "MatRetrieveValues_SeqAIJ" 2102 int MatRetrieveValues_SeqAIJ(Mat mat) 2103 { 2104 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2105 int nz = aij->i[aij->m]+aij->indexshift,ierr; 2106 2107 PetscFunctionBegin; 2108 if (aij->nonew != 1) { 2109 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2110 } 2111 if (!aij->saved_values) { 2112 SETERRQ(1,1,"Must call MatStoreValues(A);first"); 2113 } 2114 2115 /* copy values over */ 2116 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(Scalar));CHKERRQ(ierr); 2117 PetscFunctionReturn(0); 2118 } 2119 EXTERN_C_END 2120 2121 #undef __FUNC__ 2122 #define __FUNC__ "MatRetrieveValues" 2123 /*@ 2124 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2125 example, reuse of the linear part of a Jacobian, while recomputing the 2126 nonlinear portion. 2127 2128 Collect on Mat 2129 2130 Input Parameters: 2131 . mat - the matrix (currently on AIJ matrices support this option) 2132 2133 Level: advanced 2134 2135 .seealso: MatStoreValues() 2136 2137 @*/ 2138 int MatRetrieveValues(Mat mat) 2139 { 2140 int ierr,(*f)(Mat); 2141 2142 PetscFunctionBegin; 2143 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2144 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2145 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2146 2147 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void **)&f);CHKERRQ(ierr); 2148 if (f) { 2149 ierr = (*f)(mat);CHKERRQ(ierr); 2150 } else { 2151 SETERRQ(1,1,"Wrong type of matrix to retrieve values"); 2152 } 2153 PetscFunctionReturn(0); 2154 } 2155 2156 /* --------------------------------------------------------------------------------*/ 2157 2158 #undef __FUNC__ 2159 #define __FUNC__ "MatCreateSeqAIJ" 2160 /*@C 2161 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2162 (the default parallel PETSc format). For good matrix assembly performance 2163 the user should preallocate the matrix storage by setting the parameter nz 2164 (or the array nnz). By setting these parameters accurately, performance 2165 during matrix assembly can be increased by more than a factor of 50. 2166 2167 Collective on MPI_Comm 2168 2169 Input Parameters: 2170 + comm - MPI communicator, set to PETSC_COMM_SELF 2171 . m - number of rows 2172 . n - number of columns 2173 . nz - number of nonzeros per row (same for all rows) 2174 - nnz - array containing the number of nonzeros in the various rows 2175 (possibly different for each row) or PETSC_NULL 2176 2177 Output Parameter: 2178 . A - the matrix 2179 2180 Notes: 2181 The AIJ format (also called the Yale sparse matrix format or 2182 compressed row storage), is fully compatible with standard Fortran 77 2183 storage. That is, the stored row and column indices can begin at 2184 either one (as in Fortran) or zero. See the users' manual for details. 2185 2186 Specify the preallocated storage with either nz or nnz (not both). 2187 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2188 allocation. For large problems you MUST preallocate memory or you 2189 will get TERRIBLE performance, see the users' manual chapter on matrices. 2190 2191 By default, this format uses inodes (identical nodes) when possible, to 2192 improve numerical efficiency of matrix-vector products and solves. We 2193 search for consecutive rows with the same nonzero structure, thereby 2194 reusing matrix information to achieve increased efficiency. 2195 2196 Options Database Keys: 2197 + -mat_aij_no_inode - Do not use inodes 2198 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2199 - -mat_aij_oneindex - Internally use indexing starting at 1 2200 rather than 0. Note that when calling MatSetValues(), 2201 the user still MUST index entries starting at 0! 2202 2203 Level: intermediate 2204 2205 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2206 2207 @*/ 2208 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz,Mat *A) 2209 { 2210 Mat B; 2211 Mat_SeqAIJ *b; 2212 int i,len,ierr,size; 2213 PetscTruth flg; 2214 2215 PetscFunctionBegin; 2216 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2217 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Comm must be of size 1"); 2218 2219 if (nz < -2) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"nz cannot be less than -2: value %d",nz); 2220 if (nnz) { 2221 for (i=0; i<m; i++) { 2222 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,0,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2223 } 2224 } 2225 2226 *A = 0; 2227 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,"Mat",comm,MatDestroy,MatView); 2228 PLogObjectCreate(B); 2229 B->data = (void*)(b = PetscNew(Mat_SeqAIJ));CHKPTRQ(b); 2230 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2231 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2232 B->ops->destroy = MatDestroy_SeqAIJ; 2233 B->ops->view = MatView_SeqAIJ; 2234 B->factor = 0; 2235 B->lupivotthreshold = 1.0; 2236 B->mapping = 0; 2237 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2238 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2239 b->row = 0; 2240 b->col = 0; 2241 b->icol = 0; 2242 b->indexshift = 0; 2243 b->reallocs = 0; 2244 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex",&flg);CHKERRQ(ierr); 2245 if (flg) b->indexshift = -1; 2246 2247 b->m = m; B->m = m; B->M = m; 2248 b->n = n; B->n = n; B->N = n; 2249 2250 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 2251 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 2252 2253 b->imax = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(b->imax); 2254 if (!nnz) { 2255 if (nz == PETSC_DEFAULT) nz = 10; 2256 else if (nz <= 0) nz = 1; 2257 for (i=0; i<m; i++) b->imax[i] = nz; 2258 nz = nz*m; 2259 } else { 2260 nz = 0; 2261 for (i=0; i<m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2262 } 2263 2264 /* allocate the matrix space */ 2265 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 2266 b->a = (Scalar*)PetscMalloc(len);CHKPTRQ(b->a); 2267 b->j = (int*)(b->a + nz); 2268 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2269 b->i = b->j + nz; 2270 b->singlemalloc = PETSC_TRUE; 2271 b->freedata = PETSC_TRUE; 2272 2273 b->i[0] = -b->indexshift; 2274 for (i=1; i<m+1; i++) { 2275 b->i[i] = b->i[i-1] + b->imax[i-1]; 2276 } 2277 2278 /* b->ilen will count nonzeros in each row so far. */ 2279 b->ilen = (int*)PetscMalloc((m+1)*sizeof(int)); 2280 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2281 for (i=0; i<b->m; i++) { b->ilen[i] = 0;} 2282 2283 b->nz = 0; 2284 b->maxnz = nz; 2285 b->sorted = PETSC_FALSE; 2286 b->ignorezeroentries = PETSC_FALSE; 2287 b->roworiented = PETSC_TRUE; 2288 b->nonew = 0; 2289 b->diag = 0; 2290 b->solve_work = 0; 2291 b->spptr = 0; 2292 b->inode.node_count = 0; 2293 b->inode.size = 0; 2294 b->inode.limit = 5; 2295 b->inode.max_limit = 5; 2296 b->saved_values = 0; 2297 B->info.nz_unneeded = (double)b->maxnz; 2298 b->idiag = 0; 2299 b->ssor = 0; 2300 b->keepzeroedrows = PETSC_FALSE; 2301 2302 *A = B; 2303 2304 /* SuperLU is not currently supported through PETSc */ 2305 #if defined(PETSC_HAVE_SUPERLU) 2306 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu",&flg);CHKERRQ(ierr); 2307 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B);CHKERRQ(ierr); } 2308 #endif 2309 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl",&flg);CHKERRQ(ierr); 2310 if (flg) { ierr = MatUseEssl_SeqAIJ(B);CHKERRQ(ierr); } 2311 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml",&flg);CHKERRQ(ierr); 2312 if (flg) { 2313 if (!b->indexshift) SETERRQ(PETSC_ERR_LIB,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 2314 ierr = MatUseDXML_SeqAIJ(B);CHKERRQ(ierr); 2315 } 2316 ierr = OptionsHasName(PETSC_NULL,"-help",&flg);CHKERRQ(ierr); 2317 if (flg) {ierr = MatPrintHelp(B);CHKERRQ(ierr); } 2318 2319 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2320 "MatSeqAIJSetColumnIndices_SeqAIJ", 2321 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2322 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2323 "MatStoreValues_SeqAIJ", 2324 (void*)MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2325 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2326 "MatRetrieveValues_SeqAIJ", 2327 (void*)MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2328 PetscFunctionReturn(0); 2329 } 2330 2331 #undef __FUNC__ 2332 #define __FUNC__ "MatDuplicate_SeqAIJ" 2333 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2334 { 2335 Mat C; 2336 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2337 int i,len,m = a->m,shift = a->indexshift,ierr; 2338 2339 PetscFunctionBegin; 2340 *B = 0; 2341 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,"Mat",A->comm,MatDestroy,MatView); 2342 PLogObjectCreate(C); 2343 C->data = (void*)(c = PetscNew(Mat_SeqAIJ));CHKPTRQ(c); 2344 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2345 C->ops->destroy = MatDestroy_SeqAIJ; 2346 C->ops->view = MatView_SeqAIJ; 2347 C->factor = A->factor; 2348 c->row = 0; 2349 c->col = 0; 2350 c->icol = 0; 2351 c->indexshift = shift; 2352 c->keepzeroedrows = a->keepzeroedrows; 2353 C->assembled = PETSC_TRUE; 2354 2355 c->m = C->m = a->m; 2356 c->n = C->n = a->n; 2357 C->M = a->m; 2358 C->N = a->n; 2359 2360 c->imax = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(c->imax); 2361 c->ilen = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(c->ilen); 2362 for (i=0; i<m; i++) { 2363 c->imax[i] = a->imax[i]; 2364 c->ilen[i] = a->ilen[i]; 2365 } 2366 2367 /* allocate the matrix space */ 2368 c->singlemalloc = PETSC_TRUE; 2369 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 2370 c->a = (Scalar*)PetscMalloc(len);CHKPTRQ(c->a); 2371 c->j = (int*)(c->a + a->i[m] + shift); 2372 c->i = c->j + a->i[m] + shift; 2373 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2374 if (m > 0) { 2375 ierr = PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int));CHKERRQ(ierr); 2376 if (cpvalues == MAT_COPY_VALUES) { 2377 ierr = PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar));CHKERRQ(ierr); 2378 } else { 2379 ierr = PetscMemzero(c->a,(a->i[m]+shift)*sizeof(Scalar));CHKERRQ(ierr); 2380 } 2381 } 2382 2383 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2384 c->sorted = a->sorted; 2385 c->roworiented = a->roworiented; 2386 c->nonew = a->nonew; 2387 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2388 c->saved_values = 0; 2389 c->idiag = 0; 2390 c->ssor = 0; 2391 c->ignorezeroentries = a->ignorezeroentries; 2392 c->freedata = PETSC_TRUE; 2393 2394 if (a->diag) { 2395 c->diag = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(c->diag); 2396 PLogObjectMemory(C,(m+1)*sizeof(int)); 2397 for (i=0; i<m; i++) { 2398 c->diag[i] = a->diag[i]; 2399 } 2400 } else c->diag = 0; 2401 c->inode.limit = a->inode.limit; 2402 c->inode.max_limit = a->inode.max_limit; 2403 if (a->inode.size){ 2404 c->inode.size = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(c->inode.size); 2405 c->inode.node_count = a->inode.node_count; 2406 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2407 } else { 2408 c->inode.size = 0; 2409 c->inode.node_count = 0; 2410 } 2411 c->nz = a->nz; 2412 c->maxnz = a->maxnz; 2413 c->solve_work = 0; 2414 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2415 2416 *B = C; 2417 ierr = FListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2418 PetscFunctionReturn(0); 2419 } 2420 2421 #undef __FUNC__ 2422 #define __FUNC__ "MatLoad_SeqAIJ" 2423 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 2424 { 2425 Mat_SeqAIJ *a; 2426 Mat B; 2427 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N,shift; 2428 MPI_Comm comm; 2429 2430 PetscFunctionBegin; 2431 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2432 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2433 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,0,"view must have one processor"); 2434 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2435 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2436 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not matrix object in file"); 2437 M = header[1]; N = header[2]; nz = header[3]; 2438 2439 if (nz < 0) { 2440 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2441 } 2442 2443 /* read in row lengths */ 2444 rowlengths = (int*)PetscMalloc(M*sizeof(int));CHKPTRQ(rowlengths); 2445 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2446 2447 /* create our matrix */ 2448 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A);CHKERRQ(ierr); 2449 B = *A; 2450 a = (Mat_SeqAIJ*)B->data; 2451 shift = a->indexshift; 2452 2453 /* read in column indices and adjust for Fortran indexing*/ 2454 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2455 if (shift) { 2456 for (i=0; i<nz; i++) { 2457 a->j[i] += 1; 2458 } 2459 } 2460 2461 /* read in nonzero values */ 2462 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2463 2464 /* set matrix "i" values */ 2465 a->i[0] = -shift; 2466 for (i=1; i<= M; i++) { 2467 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2468 a->ilen[i-1] = rowlengths[i-1]; 2469 } 2470 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2471 2472 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2473 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2474 PetscFunctionReturn(0); 2475 } 2476 2477 #undef __FUNC__ 2478 #define __FUNC__ "MatEqual_SeqAIJ" 2479 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2480 { 2481 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2482 int ierr; 2483 2484 PetscFunctionBegin; 2485 if (B->type !=MATSEQAIJ)SETERRQ(PETSC_ERR_ARG_INCOMP,0,"Matrices must be same type"); 2486 2487 /* If the matrix dimensions are not equal,or no of nonzeros or shift */ 2488 if ((a->m != b->m) || (a->n !=b->n) ||(a->nz != b->nz)|| 2489 (a->indexshift != b->indexshift)) { 2490 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2491 } 2492 2493 /* if the a->i are the same */ 2494 ierr = PetscMemcmp(a->i,b->i,(a->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2495 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2496 2497 /* if a->j are the same */ 2498 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2499 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2500 2501 /* if a->a are the same */ 2502 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(Scalar),flg);CHKERRQ(ierr); 2503 2504 PetscFunctionReturn(0); 2505 2506 } 2507 2508 #undef __FUNC__ 2509 #define __FUNC__ "MatCreateSeqAIJWithArrays" 2510 /*@C 2511 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2512 provided by the user. 2513 2514 Coolective on MPI_Comm 2515 2516 Input Parameters: 2517 + comm - must be an MPI communicator of size 1 2518 . m - number of rows 2519 . n - number of columns 2520 . i - row indices 2521 . j - column indices 2522 - a - matrix values 2523 2524 Output Parameter: 2525 . mat - the matrix 2526 2527 Level: intermediate 2528 2529 Notes: 2530 The i, j, and a arrays are not copied by this routine, the user must free these routines 2531 once the matrix is destroyed 2532 2533 You cannot set new nonzero locations into this matrix, that will generate an error. 2534 2535 The i and j indices can be either 0- or 1 based 2536 2537 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2538 2539 @*/ 2540 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,Scalar *a,Mat *mat) 2541 { 2542 int ierr,ii; 2543 Mat_SeqAIJ *aij; 2544 2545 PetscFunctionBegin; 2546 ierr = MatCreateSeqAIJ(comm,m,n,0,0,mat);CHKERRQ(ierr); 2547 aij = (Mat_SeqAIJ*)(*mat)->data; 2548 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2549 2550 if (i[0] == 1) { 2551 aij->indexshift = -1; 2552 } else if (i[0]) { 2553 SETERRQ(1,1,"i (row indices) do not start with 0 or 1"); 2554 } 2555 aij->i = i; 2556 aij->j = j; 2557 aij->a = a; 2558 aij->singlemalloc = PETSC_FALSE; 2559 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2560 aij->freedata = PETSC_FALSE; 2561 2562 for (ii=0; ii<m; ii++) { 2563 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2564 #if defined(PETSC_BOPT_g) 2565 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]); 2566 #endif 2567 } 2568 #if defined(PETSC_BOPT_g) 2569 for (ii=0; ii<aij->i[m]; ii++) { 2570 if (j[ii] < -aij->indexshift) SETERRQ2(1,1,"Negative column index at location = %d index = %d",ii,j[ii]); 2571 if (j[ii] > n - 1 -aij->indexshift) SETERRQ2(1,1,"Column index to large at location = %d index = %d",ii,j[ii]); 2572 } 2573 #endif 2574 2575 ierr = MatAssemblyEnd_SeqAIJ(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2576 PetscFunctionReturn(0); 2577 } 2578 2579 2580 2581