1 2 /* 3 Defines the basic matrix operations for the AIJ (compressed row) 4 matrix storage format. 5 */ 6 7 #include "src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 8 #include "src/inline/spops.h" 9 #include "src/inline/dot.h" 10 #include "petscbt.h" 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 14 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *m,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 15 { 16 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 17 PetscErrorCode ierr; 18 PetscInt i,ishift; 19 20 PetscFunctionBegin; 21 *m = A->m; 22 if (!ia) PetscFunctionReturn(0); 23 ishift = 0; 24 if (symmetric && !A->structurally_symmetric) { 25 ierr = MatToSymmetricIJ_SeqAIJ(A->m,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 26 } else if (oshift == 1) { 27 PetscInt nz = a->i[A->m]; 28 /* malloc space and add 1 to i and j indices */ 29 ierr = PetscMalloc((A->m+1)*sizeof(PetscInt),ia);CHKERRQ(ierr); 30 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),ja);CHKERRQ(ierr); 31 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 32 for (i=0; i<A->m+1; i++) (*ia)[i] = a->i[i] + 1; 33 } else { 34 *ia = a->i; *ja = a->j; 35 } 36 PetscFunctionReturn(0); 37 } 38 39 #undef __FUNCT__ 40 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 41 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 42 { 43 PetscErrorCode ierr; 44 45 PetscFunctionBegin; 46 if (!ia) PetscFunctionReturn(0); 47 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 48 ierr = PetscFree(*ia);CHKERRQ(ierr); 49 ierr = PetscFree(*ja);CHKERRQ(ierr); 50 } 51 PetscFunctionReturn(0); 52 } 53 54 #undef __FUNCT__ 55 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 56 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 57 { 58 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 59 PetscErrorCode ierr; 60 PetscInt i,*collengths,*cia,*cja,n = A->n,m = A->m; 61 PetscInt nz = a->i[m],row,*jj,mr,col; 62 63 PetscFunctionBegin; 64 *nn = A->n; 65 if (!ia) PetscFunctionReturn(0); 66 if (symmetric) { 67 ierr = MatToSymmetricIJ_SeqAIJ(A->m,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 68 } else { 69 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 70 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 71 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 72 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 73 jj = a->j; 74 for (i=0; i<nz; i++) { 75 collengths[jj[i]]++; 76 } 77 cia[0] = oshift; 78 for (i=0; i<n; i++) { 79 cia[i+1] = cia[i] + collengths[i]; 80 } 81 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 82 jj = a->j; 83 for (row=0; row<m; row++) { 84 mr = a->i[row+1] - a->i[row]; 85 for (i=0; i<mr; i++) { 86 col = *jj++; 87 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 88 } 89 } 90 ierr = PetscFree(collengths);CHKERRQ(ierr); 91 *ia = cia; *ja = cja; 92 } 93 PetscFunctionReturn(0); 94 } 95 96 #undef __FUNCT__ 97 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 98 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 99 { 100 PetscErrorCode ierr; 101 102 PetscFunctionBegin; 103 if (!ia) PetscFunctionReturn(0); 104 105 ierr = PetscFree(*ia);CHKERRQ(ierr); 106 ierr = PetscFree(*ja);CHKERRQ(ierr); 107 108 PetscFunctionReturn(0); 109 } 110 111 #define CHUNKSIZE 15 112 113 #undef __FUNCT__ 114 #define __FUNCT__ "MatSetValues_SeqAIJ" 115 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 116 { 117 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 118 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 119 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 120 PetscErrorCode ierr; 121 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 122 PetscScalar *ap,value,*aa = a->a; 123 PetscTruth ignorezeroentries = ((a->ignorezeroentries && is == ADD_VALUES) ? PETSC_TRUE:PETSC_FALSE); 124 PetscTruth roworiented = a->roworiented; 125 126 PetscFunctionBegin; 127 for (k=0; k<m; k++) { /* loop over added rows */ 128 row = im[k]; 129 if (row < 0) continue; 130 #if defined(PETSC_USE_DEBUG) 131 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->m-1); 132 #endif 133 rp = aj + ai[row]; ap = aa + ai[row]; 134 rmax = imax[row]; nrow = ailen[row]; 135 low = 0; 136 high = nrow; 137 for (l=0; l<n; l++) { /* loop over added columns */ 138 if (in[l] < 0) continue; 139 #if defined(PETSC_USE_DEBUG) 140 if (in[l] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->n-1); 141 #endif 142 col = in[l]; 143 if (roworiented) { 144 value = v[l + k*n]; 145 } else { 146 value = v[k + l*m]; 147 } 148 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 149 150 if (col < lastcol) low = 0; else high = nrow; 151 lastcol = col; 152 while (high-low > 5) { 153 t = (low+high)/2; 154 if (rp[t] > col) high = t; 155 else low = t; 156 } 157 for (i=low; i<high; i++) { 158 if (rp[i] > col) break; 159 if (rp[i] == col) { 160 if (is == ADD_VALUES) ap[i] += value; 161 else ap[i] = value; 162 goto noinsert; 163 } 164 } 165 if (value == 0.0 && ignorezeroentries) goto noinsert; 166 if (nonew == 1) goto noinsert; 167 else if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 168 if (nrow >= rmax) { 169 /* there is no extra room in row, therefore enlarge */ 170 PetscInt new_nz = ai[A->m] + CHUNKSIZE,*new_i,*new_j; 171 size_t len; 172 PetscScalar *new_a; 173 174 if (nonew == -2) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix requiring new malloc()",row,col); 175 176 /* malloc new storage space */ 177 len = ((size_t) new_nz)*(sizeof(PetscInt)+sizeof(PetscScalar))+(A->m+1)*sizeof(PetscInt); 178 ierr = PetscMalloc(len,&new_a);CHKERRQ(ierr); 179 new_j = (PetscInt*)(new_a + new_nz); 180 new_i = new_j + new_nz; 181 182 /* copy over old data into new slots */ 183 for (ii=0; ii<row+1; ii++) {new_i[ii] = ai[ii];} 184 for (ii=row+1; ii<A->m+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 185 ierr = PetscMemcpy(new_j,aj,(ai[row]+nrow)*sizeof(PetscInt));CHKERRQ(ierr); 186 len = (((size_t) new_nz) - CHUNKSIZE - ai[row] - nrow ); 187 ierr = PetscMemcpy(new_j+ai[row]+nrow+CHUNKSIZE,aj+ai[row]+nrow,len*sizeof(PetscInt));CHKERRQ(ierr); 188 ierr = PetscMemcpy(new_a,aa,(((size_t) ai[row])+nrow)*sizeof(PetscScalar));CHKERRQ(ierr); 189 ierr = PetscMemcpy(new_a+ai[row]+nrow+CHUNKSIZE,aa+ai[row]+nrow,len*sizeof(PetscScalar));CHKERRQ(ierr); 190 /* free up old matrix storage */ 191 ierr = PetscFree(a->a);CHKERRQ(ierr); 192 if (!a->singlemalloc) { 193 ierr = PetscFree(a->i);CHKERRQ(ierr); 194 ierr = PetscFree(a->j);CHKERRQ(ierr); 195 } 196 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 197 a->singlemalloc = PETSC_TRUE; 198 199 rp = aj + ai[row]; ap = aa + ai[row] ; 200 rmax = imax[row] = imax[row] + CHUNKSIZE; 201 ierr = PetscLogObjectMemory(A,CHUNKSIZE*(sizeof(PetscInt) + sizeof(PetscScalar)));CHKERRQ(ierr); 202 a->maxnz += CHUNKSIZE; 203 a->reallocs++; 204 } 205 N = nrow++ - 1; a->nz++; 206 /* shift up all the later entries in this row */ 207 for (ii=N; ii>=i; ii--) { 208 rp[ii+1] = rp[ii]; 209 ap[ii+1] = ap[ii]; 210 } 211 rp[i] = col; 212 ap[i] = value; 213 noinsert:; 214 low = i + 1; 215 } 216 ailen[row] = nrow; 217 } 218 A->same_nonzero = PETSC_FALSE; 219 PetscFunctionReturn(0); 220 } 221 222 #undef __FUNCT__ 223 #define __FUNCT__ "MatGetValues_SeqAIJ" 224 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 225 { 226 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 227 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 228 PetscInt *ai = a->i,*ailen = a->ilen; 229 PetscScalar *ap,*aa = a->a,zero = 0.0; 230 231 PetscFunctionBegin; 232 for (k=0; k<m; k++) { /* loop over rows */ 233 row = im[k]; 234 if (row < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); 235 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->m-1); 236 rp = aj + ai[row]; ap = aa + ai[row]; 237 nrow = ailen[row]; 238 for (l=0; l<n; l++) { /* loop over columns */ 239 if (in[l] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); 240 if (in[l] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->n-1); 241 col = in[l] ; 242 high = nrow; low = 0; /* assume unsorted */ 243 while (high-low > 5) { 244 t = (low+high)/2; 245 if (rp[t] > col) high = t; 246 else low = t; 247 } 248 for (i=low; i<high; i++) { 249 if (rp[i] > col) break; 250 if (rp[i] == col) { 251 *v++ = ap[i]; 252 goto finished; 253 } 254 } 255 *v++ = zero; 256 finished:; 257 } 258 } 259 PetscFunctionReturn(0); 260 } 261 262 263 #undef __FUNCT__ 264 #define __FUNCT__ "MatView_SeqAIJ_Binary" 265 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 266 { 267 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 268 PetscErrorCode ierr; 269 PetscInt i,*col_lens; 270 int fd; 271 272 PetscFunctionBegin; 273 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 274 ierr = PetscMalloc((4+A->m)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 275 col_lens[0] = MAT_FILE_COOKIE; 276 col_lens[1] = A->m; 277 col_lens[2] = A->n; 278 col_lens[3] = a->nz; 279 280 /* store lengths of each row and write (including header) to file */ 281 for (i=0; i<A->m; i++) { 282 col_lens[4+i] = a->i[i+1] - a->i[i]; 283 } 284 ierr = PetscBinaryWrite(fd,col_lens,4+A->m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 285 ierr = PetscFree(col_lens);CHKERRQ(ierr); 286 287 /* store column indices (zero start index) */ 288 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 289 290 /* store nonzero values */ 291 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 292 PetscFunctionReturn(0); 293 } 294 295 EXTERN PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 296 297 #undef __FUNCT__ 298 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 299 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 300 { 301 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 302 PetscErrorCode ierr; 303 PetscInt i,j,m = A->m,shift=0; 304 char *name; 305 PetscViewerFormat format; 306 307 PetscFunctionBegin; 308 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 309 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 310 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL || format == PETSC_VIEWER_ASCII_INFO) { 311 if (a->inode.size) { 312 ierr = PetscViewerASCIIPrintf(viewer,"using I-node routines: found %D nodes, limit used is %D\n",a->inode.node_count,a->inode.limit);CHKERRQ(ierr); 313 } else { 314 ierr = PetscViewerASCIIPrintf(viewer,"not using I-node routines\n");CHKERRQ(ierr); 315 } 316 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 317 PetscInt nofinalvalue = 0; 318 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->n-!shift)) { 319 nofinalvalue = 1; 320 } 321 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 322 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->n);CHKERRQ(ierr); 323 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 324 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 325 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 326 327 for (i=0; i<m; i++) { 328 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 329 #if defined(PETSC_USE_COMPLEX) 330 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e + %18.16ei \n",i+1,a->j[j]+!shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 331 #else 332 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 333 #endif 334 } 335 } 336 if (nofinalvalue) { 337 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->n,0.0);CHKERRQ(ierr); 338 } 339 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 340 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 341 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 342 PetscFunctionReturn(0); 343 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 344 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 345 for (i=0; i<m; i++) { 346 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 347 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 348 #if defined(PETSC_USE_COMPLEX) 349 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 350 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 351 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 352 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 353 } else if (PetscRealPart(a->a[j]) != 0.0) { 354 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 355 } 356 #else 357 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 358 #endif 359 } 360 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 361 } 362 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 363 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 364 PetscInt nzd=0,fshift=1,*sptr; 365 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 366 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&sptr);CHKERRQ(ierr); 367 for (i=0; i<m; i++) { 368 sptr[i] = nzd+1; 369 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 370 if (a->j[j] >= i) { 371 #if defined(PETSC_USE_COMPLEX) 372 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 373 #else 374 if (a->a[j] != 0.0) nzd++; 375 #endif 376 } 377 } 378 } 379 sptr[m] = nzd+1; 380 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 381 for (i=0; i<m+1; i+=6) { 382 if (i+4<m) {ierr = PetscViewerASCIIPrintf(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);} 383 else if (i+3<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4]);CHKERRQ(ierr);} 384 else if (i+2<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]);CHKERRQ(ierr);} 385 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 386 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 387 else {ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr);} 388 } 389 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 390 ierr = PetscFree(sptr);CHKERRQ(ierr); 391 for (i=0; i<m; i++) { 392 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 393 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 394 } 395 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 396 } 397 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 398 for (i=0; i<m; i++) { 399 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 400 if (a->j[j] >= i) { 401 #if defined(PETSC_USE_COMPLEX) 402 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 403 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 404 } 405 #else 406 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 407 #endif 408 } 409 } 410 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 411 } 412 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 413 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 414 PetscInt cnt = 0,jcnt; 415 PetscScalar value; 416 417 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 418 for (i=0; i<m; i++) { 419 jcnt = 0; 420 for (j=0; j<A->n; j++) { 421 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 422 value = a->a[cnt++]; 423 jcnt++; 424 } else { 425 value = 0.0; 426 } 427 #if defined(PETSC_USE_COMPLEX) 428 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 429 #else 430 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 431 #endif 432 } 433 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 434 } 435 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 436 } else { 437 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 438 for (i=0; i<m; i++) { 439 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 440 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 441 #if defined(PETSC_USE_COMPLEX) 442 if (PetscImaginaryPart(a->a[j]) > 0.0) { 443 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 444 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 445 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 446 } else { 447 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 448 } 449 #else 450 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 451 #endif 452 } 453 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 454 } 455 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 456 } 457 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 458 PetscFunctionReturn(0); 459 } 460 461 #undef __FUNCT__ 462 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 463 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 464 { 465 Mat A = (Mat) Aa; 466 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 467 PetscErrorCode ierr; 468 PetscInt i,j,m = A->m,color; 469 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 470 PetscViewer viewer; 471 PetscViewerFormat format; 472 473 PetscFunctionBegin; 474 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 475 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 476 477 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 478 /* loop over matrix elements drawing boxes */ 479 480 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 481 /* Blue for negative, Cyan for zero and Red for positive */ 482 color = PETSC_DRAW_BLUE; 483 for (i=0; i<m; i++) { 484 y_l = m - i - 1.0; y_r = y_l + 1.0; 485 for (j=a->i[i]; j<a->i[i+1]; j++) { 486 x_l = a->j[j] ; x_r = x_l + 1.0; 487 #if defined(PETSC_USE_COMPLEX) 488 if (PetscRealPart(a->a[j]) >= 0.) continue; 489 #else 490 if (a->a[j] >= 0.) continue; 491 #endif 492 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 493 } 494 } 495 color = PETSC_DRAW_CYAN; 496 for (i=0; i<m; i++) { 497 y_l = m - i - 1.0; y_r = y_l + 1.0; 498 for (j=a->i[i]; j<a->i[i+1]; j++) { 499 x_l = a->j[j]; x_r = x_l + 1.0; 500 if (a->a[j] != 0.) continue; 501 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 502 } 503 } 504 color = PETSC_DRAW_RED; 505 for (i=0; i<m; i++) { 506 y_l = m - i - 1.0; y_r = y_l + 1.0; 507 for (j=a->i[i]; j<a->i[i+1]; j++) { 508 x_l = a->j[j]; x_r = x_l + 1.0; 509 #if defined(PETSC_USE_COMPLEX) 510 if (PetscRealPart(a->a[j]) <= 0.) continue; 511 #else 512 if (a->a[j] <= 0.) continue; 513 #endif 514 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 515 } 516 } 517 } else { 518 /* use contour shading to indicate magnitude of values */ 519 /* first determine max of all nonzero values */ 520 PetscInt nz = a->nz,count; 521 PetscDraw popup; 522 PetscReal scale; 523 524 for (i=0; i<nz; i++) { 525 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 526 } 527 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 528 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 529 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 530 count = 0; 531 for (i=0; i<m; i++) { 532 y_l = m - i - 1.0; y_r = y_l + 1.0; 533 for (j=a->i[i]; j<a->i[i+1]; j++) { 534 x_l = a->j[j]; x_r = x_l + 1.0; 535 color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count])); 536 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 537 count++; 538 } 539 } 540 } 541 PetscFunctionReturn(0); 542 } 543 544 #undef __FUNCT__ 545 #define __FUNCT__ "MatView_SeqAIJ_Draw" 546 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 547 { 548 PetscErrorCode ierr; 549 PetscDraw draw; 550 PetscReal xr,yr,xl,yl,h,w; 551 PetscTruth isnull; 552 553 PetscFunctionBegin; 554 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 555 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 556 if (isnull) PetscFunctionReturn(0); 557 558 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 559 xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0; 560 xr += w; yr += h; xl = -w; yl = -h; 561 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 562 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 563 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 564 PetscFunctionReturn(0); 565 } 566 567 #undef __FUNCT__ 568 #define __FUNCT__ "MatView_SeqAIJ" 569 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 570 { 571 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 572 PetscErrorCode ierr; 573 PetscTruth issocket,iascii,isbinary,isdraw; 574 575 PetscFunctionBegin; 576 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 577 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 578 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 579 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 580 if (issocket) { 581 ierr = PetscViewerSocketPutSparse_Private(viewer,A->m,A->n,a->nz,a->a,a->i,a->j);CHKERRQ(ierr); 582 } else if (iascii) { 583 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 584 } else if (isbinary) { 585 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 586 } else if (isdraw) { 587 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 588 } else { 589 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 590 } 591 PetscFunctionReturn(0); 592 } 593 594 #undef __FUNCT__ 595 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 596 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 597 { 598 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 599 PetscErrorCode ierr; 600 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 601 PetscInt m = A->m,*ip,N,*ailen = a->ilen,rmax = 0; 602 PetscScalar *aa = a->a,*ap; 603 PetscReal ratio=0.6; 604 605 PetscFunctionBegin; 606 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 607 608 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 609 for (i=1; i<m; i++) { 610 /* move each row back by the amount of empty slots (fshift) before it*/ 611 fshift += imax[i-1] - ailen[i-1]; 612 rmax = PetscMax(rmax,ailen[i]); 613 if (fshift) { 614 ip = aj + ai[i] ; 615 ap = aa + ai[i] ; 616 N = ailen[i]; 617 for (j=0; j<N; j++) { 618 ip[j-fshift] = ip[j]; 619 ap[j-fshift] = ap[j]; 620 } 621 } 622 ai[i] = ai[i-1] + ailen[i-1]; 623 } 624 if (m) { 625 fshift += imax[m-1] - ailen[m-1]; 626 ai[m] = ai[m-1] + ailen[m-1]; 627 } 628 /* reset ilen and imax for each row */ 629 for (i=0; i<m; i++) { 630 ailen[i] = imax[i] = ai[i+1] - ai[i]; 631 } 632 a->nz = ai[m]; 633 634 /* diagonals may have moved, so kill the diagonal pointers */ 635 if (fshift && a->diag) { 636 ierr = PetscFree(a->diag);CHKERRQ(ierr); 637 ierr = PetscLogObjectMemory(A,-(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 638 a->diag = 0; 639 } 640 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->n,fshift,a->nz); 641 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues() is %D\n",a->reallocs); 642 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Maximum nonzeros in any row is %D\n",rmax); 643 a->reallocs = 0; 644 A->info.nz_unneeded = (double)fshift; 645 a->rmax = rmax; 646 647 /* check for identical nodes. If found, use inode functions */ 648 if (!a->inode.checked){ 649 ierr = Mat_AIJ_CheckInode(A,(PetscTruth)(!fshift));CHKERRQ(ierr); 650 } 651 652 /* check for zero rows. If found a large number of zero rows, use CompressedRow functions */ 653 if (!a->inode.use && a->compressedrow.use){ 654 ierr = Mat_CheckCompressedRow(A,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 655 } 656 A->same_nonzero = PETSC_TRUE; 657 PetscFunctionReturn(0); 658 } 659 660 #undef __FUNCT__ 661 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 662 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 663 { 664 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 665 PetscErrorCode ierr; 666 667 PetscFunctionBegin; 668 ierr = PetscMemzero(a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 669 PetscFunctionReturn(0); 670 } 671 672 #undef __FUNCT__ 673 #define __FUNCT__ "MatDestroy_SeqAIJ" 674 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 675 { 676 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 677 PetscErrorCode ierr; 678 679 PetscFunctionBegin; 680 #if defined(PETSC_USE_LOG) 681 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->m,A->n,a->nz); 682 #endif 683 if (a->freedata) { 684 ierr = PetscFree(a->a);CHKERRQ(ierr); 685 if (!a->singlemalloc) { 686 ierr = PetscFree(a->i);CHKERRQ(ierr); 687 ierr = PetscFree(a->j);CHKERRQ(ierr); 688 } 689 } 690 if (a->row) { 691 ierr = ISDestroy(a->row);CHKERRQ(ierr); 692 } 693 if (a->col) { 694 ierr = ISDestroy(a->col);CHKERRQ(ierr); 695 } 696 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 697 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 698 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 699 if (a->idiag) {ierr = PetscFree(a->idiag);CHKERRQ(ierr);} 700 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 701 if (a->inode.size) {ierr = PetscFree(a->inode.size);CHKERRQ(ierr);} 702 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 703 if (a->saved_values) {ierr = PetscFree(a->saved_values);CHKERRQ(ierr);} 704 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 705 if (a->xtoy) {ierr = PetscFree(a->xtoy);CHKERRQ(ierr);} 706 if (a->compressedrow.use){ierr = PetscFree(a->compressedrow.i);} 707 708 ierr = PetscFree(a);CHKERRQ(ierr); 709 710 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 711 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 712 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 713 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqsbaij_C","",PETSC_NULL);CHKERRQ(ierr); 714 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 715 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatIsTranspose_C","",PETSC_NULL);CHKERRQ(ierr); 716 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 717 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatReorderForNonzeroDiagonal_C","",PETSC_NULL);CHKERRQ(ierr); 718 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatAdjustForInodes_C","",PETSC_NULL);CHKERRQ(ierr); 719 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJGetInodeSizes_C","",PETSC_NULL);CHKERRQ(ierr); 720 PetscFunctionReturn(0); 721 } 722 723 #undef __FUNCT__ 724 #define __FUNCT__ "MatCompress_SeqAIJ" 725 PetscErrorCode MatCompress_SeqAIJ(Mat A) 726 { 727 PetscFunctionBegin; 728 PetscFunctionReturn(0); 729 } 730 731 #undef __FUNCT__ 732 #define __FUNCT__ "MatSetOption_SeqAIJ" 733 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op) 734 { 735 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 736 737 PetscFunctionBegin; 738 switch (op) { 739 case MAT_ROW_ORIENTED: 740 a->roworiented = PETSC_TRUE; 741 break; 742 case MAT_KEEP_ZEROED_ROWS: 743 a->keepzeroedrows = PETSC_TRUE; 744 break; 745 case MAT_COLUMN_ORIENTED: 746 a->roworiented = PETSC_FALSE; 747 break; 748 case MAT_COLUMNS_SORTED: 749 a->sorted = PETSC_TRUE; 750 break; 751 case MAT_COLUMNS_UNSORTED: 752 a->sorted = PETSC_FALSE; 753 break; 754 case MAT_NO_NEW_NONZERO_LOCATIONS: 755 a->nonew = 1; 756 break; 757 case MAT_NEW_NONZERO_LOCATION_ERR: 758 a->nonew = -1; 759 break; 760 case MAT_NEW_NONZERO_ALLOCATION_ERR: 761 a->nonew = -2; 762 break; 763 case MAT_YES_NEW_NONZERO_LOCATIONS: 764 a->nonew = 0; 765 break; 766 case MAT_IGNORE_ZERO_ENTRIES: 767 a->ignorezeroentries = PETSC_TRUE; 768 break; 769 case MAT_USE_INODES: 770 a->inode.use = PETSC_TRUE; 771 break; 772 case MAT_DO_NOT_USE_INODES: 773 a->inode.use = PETSC_FALSE; 774 break; 775 case MAT_USE_COMPRESSEDROW: 776 a->compressedrow.use = PETSC_TRUE; 777 break; 778 case MAT_DO_NOT_USE_COMPRESSEDROW: 779 a->compressedrow.use = PETSC_FALSE; 780 break; 781 case MAT_ROWS_SORTED: 782 case MAT_ROWS_UNSORTED: 783 case MAT_YES_NEW_DIAGONALS: 784 case MAT_IGNORE_OFF_PROC_ENTRIES: 785 case MAT_USE_HASH_TABLE: 786 PetscLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 787 break; 788 case MAT_NO_NEW_DIAGONALS: 789 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 790 case MAT_INODE_LIMIT_1: 791 a->inode.limit = 1; 792 break; 793 case MAT_INODE_LIMIT_2: 794 a->inode.limit = 2; 795 break; 796 case MAT_INODE_LIMIT_3: 797 a->inode.limit = 3; 798 break; 799 case MAT_INODE_LIMIT_4: 800 a->inode.limit = 4; 801 break; 802 case MAT_INODE_LIMIT_5: 803 a->inode.limit = 5; 804 break; 805 case MAT_SYMMETRIC: 806 case MAT_STRUCTURALLY_SYMMETRIC: 807 case MAT_NOT_SYMMETRIC: 808 case MAT_NOT_STRUCTURALLY_SYMMETRIC: 809 case MAT_HERMITIAN: 810 case MAT_NOT_HERMITIAN: 811 case MAT_SYMMETRY_ETERNAL: 812 case MAT_NOT_SYMMETRY_ETERNAL: 813 break; 814 default: 815 SETERRQ(PETSC_ERR_SUP,"unknown option"); 816 } 817 PetscFunctionReturn(0); 818 } 819 820 #undef __FUNCT__ 821 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 822 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 823 { 824 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 825 PetscErrorCode ierr; 826 PetscInt i,j,n; 827 PetscScalar *x,zero = 0.0; 828 829 PetscFunctionBegin; 830 ierr = VecSet(&zero,v);CHKERRQ(ierr); 831 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 832 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 833 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 834 for (i=0; i<A->m; i++) { 835 for (j=a->i[i]; j<a->i[i+1]; j++) { 836 if (a->j[j] == i) { 837 x[i] = a->a[j]; 838 break; 839 } 840 } 841 } 842 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 843 PetscFunctionReturn(0); 844 } 845 846 #undef __FUNCT__ 847 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 848 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 849 { 850 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 851 PetscScalar *x,*y; 852 PetscErrorCode ierr; 853 PetscInt m = A->m; 854 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 855 PetscScalar *v,alpha; 856 PetscInt n,i,*idx,*ii,*ridx=PETSC_NULL; 857 Mat_CompressedRow cprow = a->compressedrow; 858 PetscTruth usecprow = cprow.use; 859 #endif 860 861 PetscFunctionBegin; 862 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 863 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 864 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 865 866 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 867 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 868 #else 869 if (usecprow){ 870 m = cprow.nrows; 871 ii = cprow.i; 872 ridx = cprow.rindex; 873 } else { 874 ii = a->i; 875 } 876 for (i=0; i<m; i++) { 877 idx = a->j + ii[i] ; 878 v = a->a + ii[i] ; 879 n = ii[i+1] - ii[i]; 880 if (usecprow){ 881 alpha = x[ridx[i]]; 882 } else { 883 alpha = x[i]; 884 } 885 while (n-->0) {y[*idx++] += alpha * *v++;} 886 } 887 #endif 888 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 889 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 890 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 891 PetscFunctionReturn(0); 892 } 893 894 #undef __FUNCT__ 895 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 896 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 897 { 898 PetscScalar zero = 0.0; 899 PetscErrorCode ierr; 900 901 PetscFunctionBegin; 902 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 903 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 904 PetscFunctionReturn(0); 905 } 906 907 908 #undef __FUNCT__ 909 #define __FUNCT__ "MatMult_SeqAIJ" 910 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 911 { 912 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 913 PetscScalar *x,*y,*aa; 914 PetscErrorCode ierr; 915 PetscInt m=A->m,*aj,*ii; 916 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 917 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 918 PetscScalar sum; 919 PetscTruth usecprow=a->compressedrow.use; 920 #endif 921 922 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 923 #pragma disjoint(*x,*y,*aa) 924 #endif 925 926 PetscFunctionBegin; 927 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 928 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 929 aj = a->j; 930 aa = a->a; 931 ii = a->i; 932 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 933 fortranmultaij_(&m,x,ii,aj,aa,y); 934 #else 935 if (usecprow){ /* use compressed row format */ 936 m = a->compressedrow.nrows; 937 ii = a->compressedrow.i; 938 ridx = a->compressedrow.rindex; 939 for (i=0; i<m; i++){ 940 n = ii[i+1] - ii[i]; 941 aj = a->j + ii[i]; 942 aa = a->a + ii[i]; 943 sum = 0.0; 944 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 945 y[*ridx++] = sum; 946 } 947 } else { /* do not use compressed row format */ 948 for (i=0; i<m; i++) { 949 jrow = ii[i]; 950 n = ii[i+1] - jrow; 951 sum = 0.0; 952 for (j=0; j<n; j++) { 953 sum += aa[jrow]*x[aj[jrow]]; jrow++; 954 } 955 y[i] = sum; 956 } 957 } 958 #endif 959 ierr = PetscLogFlops(2*a->nz - m);CHKERRQ(ierr); 960 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 961 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 962 PetscFunctionReturn(0); 963 } 964 965 #undef __FUNCT__ 966 #define __FUNCT__ "MatMultAdd_SeqAIJ" 967 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 968 { 969 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 970 PetscScalar *x,*y,*z,*aa; 971 PetscErrorCode ierr; 972 PetscInt m = A->m,*aj,*ii; 973 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 974 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 975 PetscScalar sum; 976 PetscTruth usecprow=a->compressedrow.use; 977 #endif 978 979 PetscFunctionBegin; 980 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 981 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 982 if (zz != yy) { 983 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 984 } else { 985 z = y; 986 } 987 988 aj = a->j; 989 aa = a->a; 990 ii = a->i; 991 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 992 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 993 #else 994 if (usecprow){ /* use compressed row format */ 995 if (zz != yy){ 996 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 997 } 998 m = a->compressedrow.nrows; 999 ii = a->compressedrow.i; 1000 ridx = a->compressedrow.rindex; 1001 for (i=0; i<m; i++){ 1002 n = ii[i+1] - ii[i]; 1003 aj = a->j + ii[i]; 1004 aa = a->a + ii[i]; 1005 sum = y[*ridx]; 1006 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 1007 z[*ridx++] = sum; 1008 } 1009 } else { /* do not use compressed row format */ 1010 for (i=0; i<m; i++) { 1011 jrow = ii[i]; 1012 n = ii[i+1] - jrow; 1013 sum = y[i]; 1014 for (j=0; j<n; j++) { 1015 sum += aa[jrow]*x[aj[jrow]]; jrow++; 1016 } 1017 z[i] = sum; 1018 } 1019 } 1020 #endif 1021 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 1022 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1023 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1024 if (zz != yy) { 1025 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1026 } 1027 PetscFunctionReturn(0); 1028 } 1029 1030 /* 1031 Adds diagonal pointers to sparse matrix structure. 1032 */ 1033 #undef __FUNCT__ 1034 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1035 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1036 { 1037 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1038 PetscErrorCode ierr; 1039 PetscInt i,j,*diag,m = A->m; 1040 1041 PetscFunctionBegin; 1042 if (a->diag) PetscFunctionReturn(0); 1043 1044 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&diag);CHKERRQ(ierr); 1045 ierr = PetscLogObjectMemory(A,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 1046 for (i=0; i<A->m; i++) { 1047 diag[i] = a->i[i+1]; 1048 for (j=a->i[i]; j<a->i[i+1]; j++) { 1049 if (a->j[j] == i) { 1050 diag[i] = j; 1051 break; 1052 } 1053 } 1054 } 1055 a->diag = diag; 1056 PetscFunctionReturn(0); 1057 } 1058 1059 /* 1060 Checks for missing diagonals 1061 */ 1062 #undef __FUNCT__ 1063 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1064 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A) 1065 { 1066 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1067 PetscErrorCode ierr; 1068 PetscInt *diag,*jj = a->j,i; 1069 1070 PetscFunctionBegin; 1071 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1072 diag = a->diag; 1073 for (i=0; i<A->m; i++) { 1074 if (jj[diag[i]] != i) { 1075 SETERRQ1(PETSC_ERR_PLIB,"Matrix is missing diagonal number %D",i); 1076 } 1077 } 1078 PetscFunctionReturn(0); 1079 } 1080 1081 #undef __FUNCT__ 1082 #define __FUNCT__ "MatRelax_SeqAIJ" 1083 PetscErrorCode MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1084 { 1085 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1086 PetscScalar *x,d,*xs,sum,*t,scale,*idiag=0,*mdiag; 1087 const PetscScalar *v = a->a, *b, *bs,*xb, *ts; 1088 PetscErrorCode ierr; 1089 PetscInt n = A->n,m = A->m,i; 1090 const PetscInt *idx,*diag; 1091 1092 PetscFunctionBegin; 1093 its = its*lits; 1094 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 1095 1096 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 1097 diag = a->diag; 1098 if (!a->idiag) { 1099 ierr = PetscMalloc(3*m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 1100 a->ssor = a->idiag + m; 1101 mdiag = a->ssor + m; 1102 1103 v = a->a; 1104 1105 /* this is wrong when fshift omega changes each iteration */ 1106 if (omega == 1.0 && !fshift) { 1107 for (i=0; i<m; i++) { 1108 mdiag[i] = v[diag[i]]; 1109 a->idiag[i] = 1.0/v[diag[i]]; 1110 } 1111 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1112 } else { 1113 for (i=0; i<m; i++) { 1114 mdiag[i] = v[diag[i]]; 1115 a->idiag[i] = omega/(fshift + v[diag[i]]); 1116 } 1117 ierr = PetscLogFlops(2*m);CHKERRQ(ierr); 1118 } 1119 } 1120 t = a->ssor; 1121 idiag = a->idiag; 1122 mdiag = a->idiag + 2*m; 1123 1124 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1125 if (xx != bb) { 1126 ierr = VecGetArray(bb,(PetscScalar**)&b);CHKERRQ(ierr); 1127 } else { 1128 b = x; 1129 } 1130 1131 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1132 xs = x; 1133 if (flag == SOR_APPLY_UPPER) { 1134 /* apply (U + D/omega) to the vector */ 1135 bs = b; 1136 for (i=0; i<m; i++) { 1137 d = fshift + a->a[diag[i]]; 1138 n = a->i[i+1] - diag[i] - 1; 1139 idx = a->j + diag[i] + 1; 1140 v = a->a + diag[i] + 1; 1141 sum = b[i]*d/omega; 1142 SPARSEDENSEDOT(sum,bs,v,idx,n); 1143 x[i] = sum; 1144 } 1145 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1146 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1147 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1148 PetscFunctionReturn(0); 1149 } 1150 1151 1152 /* Let A = L + U + D; where L is lower trianglar, 1153 U is upper triangular, E is diagonal; This routine applies 1154 1155 (L + E)^{-1} A (U + E)^{-1} 1156 1157 to a vector efficiently using Eisenstat's trick. This is for 1158 the case of SSOR preconditioner, so E is D/omega where omega 1159 is the relaxation factor. 1160 */ 1161 1162 if (flag == SOR_APPLY_LOWER) { 1163 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1164 } else if (flag & SOR_EISENSTAT) { 1165 /* Let A = L + U + D; where L is lower trianglar, 1166 U is upper triangular, E is diagonal; This routine applies 1167 1168 (L + E)^{-1} A (U + E)^{-1} 1169 1170 to a vector efficiently using Eisenstat's trick. This is for 1171 the case of SSOR preconditioner, so E is D/omega where omega 1172 is the relaxation factor. 1173 */ 1174 scale = (2.0/omega) - 1.0; 1175 1176 /* x = (E + U)^{-1} b */ 1177 for (i=m-1; i>=0; i--) { 1178 n = a->i[i+1] - diag[i] - 1; 1179 idx = a->j + diag[i] + 1; 1180 v = a->a + diag[i] + 1; 1181 sum = b[i]; 1182 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1183 x[i] = sum*idiag[i]; 1184 } 1185 1186 /* t = b - (2*E - D)x */ 1187 v = a->a; 1188 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1189 1190 /* t = (E + L)^{-1}t */ 1191 ts = t; 1192 diag = a->diag; 1193 for (i=0; i<m; i++) { 1194 n = diag[i] - a->i[i]; 1195 idx = a->j + a->i[i]; 1196 v = a->a + a->i[i]; 1197 sum = t[i]; 1198 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1199 t[i] = sum*idiag[i]; 1200 /* x = x + t */ 1201 x[i] += t[i]; 1202 } 1203 1204 ierr = PetscLogFlops(6*m-1 + 2*a->nz);CHKERRQ(ierr); 1205 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1206 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1207 PetscFunctionReturn(0); 1208 } 1209 if (flag & SOR_ZERO_INITIAL_GUESS) { 1210 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1211 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1212 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)b); 1213 #else 1214 for (i=0; i<m; i++) { 1215 n = diag[i] - a->i[i]; 1216 idx = a->j + a->i[i]; 1217 v = a->a + a->i[i]; 1218 sum = b[i]; 1219 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1220 x[i] = sum*idiag[i]; 1221 } 1222 #endif 1223 xb = x; 1224 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1225 } else xb = b; 1226 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1227 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1228 for (i=0; i<m; i++) { 1229 x[i] *= mdiag[i]; 1230 } 1231 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1232 } 1233 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1234 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1235 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)xb); 1236 #else 1237 for (i=m-1; i>=0; i--) { 1238 n = a->i[i+1] - diag[i] - 1; 1239 idx = a->j + diag[i] + 1; 1240 v = a->a + diag[i] + 1; 1241 sum = xb[i]; 1242 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1243 x[i] = sum*idiag[i]; 1244 } 1245 #endif 1246 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1247 } 1248 its--; 1249 } 1250 while (its--) { 1251 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1252 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1253 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1254 #else 1255 for (i=0; i<m; i++) { 1256 d = fshift + a->a[diag[i]]; 1257 n = a->i[i+1] - a->i[i]; 1258 idx = a->j + a->i[i]; 1259 v = a->a + a->i[i]; 1260 sum = b[i]; 1261 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1262 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1263 } 1264 #endif 1265 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1266 } 1267 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1268 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1269 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1270 #else 1271 for (i=m-1; i>=0; i--) { 1272 d = fshift + a->a[diag[i]]; 1273 n = a->i[i+1] - a->i[i]; 1274 idx = a->j + a->i[i]; 1275 v = a->a + a->i[i]; 1276 sum = b[i]; 1277 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1278 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1279 } 1280 #endif 1281 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1282 } 1283 } 1284 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1285 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1286 PetscFunctionReturn(0); 1287 } 1288 1289 #undef __FUNCT__ 1290 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1291 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1292 { 1293 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1294 1295 PetscFunctionBegin; 1296 info->rows_global = (double)A->m; 1297 info->columns_global = (double)A->n; 1298 info->rows_local = (double)A->m; 1299 info->columns_local = (double)A->n; 1300 info->block_size = 1.0; 1301 info->nz_allocated = (double)a->maxnz; 1302 info->nz_used = (double)a->nz; 1303 info->nz_unneeded = (double)(a->maxnz - a->nz); 1304 info->assemblies = (double)A->num_ass; 1305 info->mallocs = (double)a->reallocs; 1306 info->memory = A->mem; 1307 if (A->factor) { 1308 info->fill_ratio_given = A->info.fill_ratio_given; 1309 info->fill_ratio_needed = A->info.fill_ratio_needed; 1310 info->factor_mallocs = A->info.factor_mallocs; 1311 } else { 1312 info->fill_ratio_given = 0; 1313 info->fill_ratio_needed = 0; 1314 info->factor_mallocs = 0; 1315 } 1316 PetscFunctionReturn(0); 1317 } 1318 1319 #undef __FUNCT__ 1320 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1321 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,IS is,const PetscScalar *diag) 1322 { 1323 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1324 PetscErrorCode ierr; 1325 PetscInt i,N,*rows,m = A->m - 1; 1326 1327 PetscFunctionBegin; 1328 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1329 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1330 if (a->keepzeroedrows) { 1331 for (i=0; i<N; i++) { 1332 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1333 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1334 } 1335 if (diag) { 1336 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1337 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1338 for (i=0; i<N; i++) { 1339 a->a[a->diag[rows[i]]] = *diag; 1340 } 1341 } 1342 A->same_nonzero = PETSC_TRUE; 1343 } else { 1344 if (diag) { 1345 for (i=0; i<N; i++) { 1346 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1347 if (a->ilen[rows[i]] > 0) { 1348 a->ilen[rows[i]] = 1; 1349 a->a[a->i[rows[i]]] = *diag; 1350 a->j[a->i[rows[i]]] = rows[i]; 1351 } else { /* in case row was completely empty */ 1352 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1353 } 1354 } 1355 } else { 1356 for (i=0; i<N; i++) { 1357 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1358 a->ilen[rows[i]] = 0; 1359 } 1360 } 1361 A->same_nonzero = PETSC_FALSE; 1362 } 1363 ierr = ISRestoreIndices(is,&rows);CHKERRQ(ierr); 1364 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1365 PetscFunctionReturn(0); 1366 } 1367 1368 #undef __FUNCT__ 1369 #define __FUNCT__ "MatGetRow_SeqAIJ" 1370 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1371 { 1372 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1373 PetscInt *itmp; 1374 1375 PetscFunctionBegin; 1376 if (row < 0 || row >= A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1377 1378 *nz = a->i[row+1] - a->i[row]; 1379 if (v) *v = a->a + a->i[row]; 1380 if (idx) { 1381 itmp = a->j + a->i[row]; 1382 if (*nz) { 1383 *idx = itmp; 1384 } 1385 else *idx = 0; 1386 } 1387 PetscFunctionReturn(0); 1388 } 1389 1390 /* remove this function? */ 1391 #undef __FUNCT__ 1392 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1393 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1394 { 1395 PetscFunctionBegin; 1396 PetscFunctionReturn(0); 1397 } 1398 1399 #undef __FUNCT__ 1400 #define __FUNCT__ "MatNorm_SeqAIJ" 1401 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1402 { 1403 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1404 PetscScalar *v = a->a; 1405 PetscReal sum = 0.0; 1406 PetscErrorCode ierr; 1407 PetscInt i,j; 1408 1409 PetscFunctionBegin; 1410 if (type == NORM_FROBENIUS) { 1411 for (i=0; i<a->nz; i++) { 1412 #if defined(PETSC_USE_COMPLEX) 1413 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1414 #else 1415 sum += (*v)*(*v); v++; 1416 #endif 1417 } 1418 *nrm = sqrt(sum); 1419 } else if (type == NORM_1) { 1420 PetscReal *tmp; 1421 PetscInt *jj = a->j; 1422 ierr = PetscMalloc((A->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1423 ierr = PetscMemzero(tmp,A->n*sizeof(PetscReal));CHKERRQ(ierr); 1424 *nrm = 0.0; 1425 for (j=0; j<a->nz; j++) { 1426 tmp[*jj++] += PetscAbsScalar(*v); v++; 1427 } 1428 for (j=0; j<A->n; j++) { 1429 if (tmp[j] > *nrm) *nrm = tmp[j]; 1430 } 1431 ierr = PetscFree(tmp);CHKERRQ(ierr); 1432 } else if (type == NORM_INFINITY) { 1433 *nrm = 0.0; 1434 for (j=0; j<A->m; j++) { 1435 v = a->a + a->i[j]; 1436 sum = 0.0; 1437 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1438 sum += PetscAbsScalar(*v); v++; 1439 } 1440 if (sum > *nrm) *nrm = sum; 1441 } 1442 } else { 1443 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1444 } 1445 PetscFunctionReturn(0); 1446 } 1447 1448 #undef __FUNCT__ 1449 #define __FUNCT__ "MatTranspose_SeqAIJ" 1450 PetscErrorCode MatTranspose_SeqAIJ(Mat A,Mat *B) 1451 { 1452 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1453 Mat C; 1454 PetscErrorCode ierr; 1455 PetscInt i,*aj = a->j,*ai = a->i,m = A->m,len,*col; 1456 PetscScalar *array = a->a; 1457 1458 PetscFunctionBegin; 1459 if (!B && m != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1460 ierr = PetscMalloc((1+A->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1461 ierr = PetscMemzero(col,(1+A->n)*sizeof(PetscInt));CHKERRQ(ierr); 1462 1463 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1464 ierr = MatCreate(A->comm,A->n,m,A->n,m,&C);CHKERRQ(ierr); 1465 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1466 ierr = MatSeqAIJSetPreallocation(C,0,col);CHKERRQ(ierr); 1467 ierr = PetscFree(col);CHKERRQ(ierr); 1468 for (i=0; i<m; i++) { 1469 len = ai[i+1]-ai[i]; 1470 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1471 array += len; 1472 aj += len; 1473 } 1474 1475 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1476 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1477 1478 if (B) { 1479 *B = C; 1480 } else { 1481 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1482 } 1483 PetscFunctionReturn(0); 1484 } 1485 1486 EXTERN_C_BEGIN 1487 #undef __FUNCT__ 1488 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1489 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1490 { 1491 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1492 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1493 PetscErrorCode ierr; 1494 PetscInt ma,na,mb,nb, i; 1495 1496 PetscFunctionBegin; 1497 bij = (Mat_SeqAIJ *) B->data; 1498 1499 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1500 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1501 if (ma!=nb || na!=mb){ 1502 *f = PETSC_FALSE; 1503 PetscFunctionReturn(0); 1504 } 1505 aii = aij->i; bii = bij->i; 1506 adx = aij->j; bdx = bij->j; 1507 va = aij->a; vb = bij->a; 1508 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1509 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1510 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1511 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1512 1513 *f = PETSC_TRUE; 1514 for (i=0; i<ma; i++) { 1515 while (aptr[i]<aii[i+1]) { 1516 PetscInt idc,idr; 1517 PetscScalar vc,vr; 1518 /* column/row index/value */ 1519 idc = adx[aptr[i]]; 1520 idr = bdx[bptr[idc]]; 1521 vc = va[aptr[i]]; 1522 vr = vb[bptr[idc]]; 1523 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1524 *f = PETSC_FALSE; 1525 goto done; 1526 } else { 1527 aptr[i]++; 1528 if (B || i!=idc) bptr[idc]++; 1529 } 1530 } 1531 } 1532 done: 1533 ierr = PetscFree(aptr);CHKERRQ(ierr); 1534 if (B) { 1535 ierr = PetscFree(bptr);CHKERRQ(ierr); 1536 } 1537 PetscFunctionReturn(0); 1538 } 1539 EXTERN_C_END 1540 1541 #undef __FUNCT__ 1542 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1543 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1544 { 1545 PetscErrorCode ierr; 1546 PetscFunctionBegin; 1547 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1548 PetscFunctionReturn(0); 1549 } 1550 1551 #undef __FUNCT__ 1552 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1553 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1554 { 1555 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1556 PetscScalar *l,*r,x,*v; 1557 PetscErrorCode ierr; 1558 PetscInt i,j,m = A->m,n = A->n,M,nz = a->nz,*jj; 1559 1560 PetscFunctionBegin; 1561 if (ll) { 1562 /* The local size is used so that VecMPI can be passed to this routine 1563 by MatDiagonalScale_MPIAIJ */ 1564 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1565 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1566 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1567 v = a->a; 1568 for (i=0; i<m; i++) { 1569 x = l[i]; 1570 M = a->i[i+1] - a->i[i]; 1571 for (j=0; j<M; j++) { (*v++) *= x;} 1572 } 1573 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1574 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1575 } 1576 if (rr) { 1577 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1578 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1579 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1580 v = a->a; jj = a->j; 1581 for (i=0; i<nz; i++) { 1582 (*v++) *= r[*jj++]; 1583 } 1584 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1585 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1586 } 1587 PetscFunctionReturn(0); 1588 } 1589 1590 #undef __FUNCT__ 1591 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1592 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1593 { 1594 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1595 PetscErrorCode ierr; 1596 PetscInt *smap,i,k,kstart,kend,oldcols = A->n,*lens; 1597 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1598 PetscInt *irow,*icol,nrows,ncols; 1599 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1600 PetscScalar *a_new,*mat_a; 1601 Mat C; 1602 PetscTruth stride; 1603 1604 PetscFunctionBegin; 1605 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1606 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1607 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1608 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1609 1610 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1611 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1612 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1613 1614 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1615 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1616 if (stride && step == 1) { 1617 /* special case of contiguous rows */ 1618 ierr = PetscMalloc((2*nrows+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1619 starts = lens + nrows; 1620 /* loop over new rows determining lens and starting points */ 1621 for (i=0; i<nrows; i++) { 1622 kstart = ai[irow[i]]; 1623 kend = kstart + ailen[irow[i]]; 1624 for (k=kstart; k<kend; k++) { 1625 if (aj[k] >= first) { 1626 starts[i] = k; 1627 break; 1628 } 1629 } 1630 sum = 0; 1631 while (k < kend) { 1632 if (aj[k++] >= first+ncols) break; 1633 sum++; 1634 } 1635 lens[i] = sum; 1636 } 1637 /* create submatrix */ 1638 if (scall == MAT_REUSE_MATRIX) { 1639 PetscInt n_cols,n_rows; 1640 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1641 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1642 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1643 C = *B; 1644 } else { 1645 ierr = MatCreate(A->comm,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE,&C);CHKERRQ(ierr); 1646 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1647 ierr = MatSeqAIJSetPreallocation(C,0,lens);CHKERRQ(ierr); 1648 } 1649 c = (Mat_SeqAIJ*)C->data; 1650 1651 /* loop over rows inserting into submatrix */ 1652 a_new = c->a; 1653 j_new = c->j; 1654 i_new = c->i; 1655 1656 for (i=0; i<nrows; i++) { 1657 ii = starts[i]; 1658 lensi = lens[i]; 1659 for (k=0; k<lensi; k++) { 1660 *j_new++ = aj[ii+k] - first; 1661 } 1662 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1663 a_new += lensi; 1664 i_new[i+1] = i_new[i] + lensi; 1665 c->ilen[i] = lensi; 1666 } 1667 ierr = PetscFree(lens);CHKERRQ(ierr); 1668 } else { 1669 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1670 ierr = PetscMalloc((1+oldcols)*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1671 1672 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1673 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1674 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1675 /* determine lens of each row */ 1676 for (i=0; i<nrows; i++) { 1677 kstart = ai[irow[i]]; 1678 kend = kstart + a->ilen[irow[i]]; 1679 lens[i] = 0; 1680 for (k=kstart; k<kend; k++) { 1681 if (smap[aj[k]]) { 1682 lens[i]++; 1683 } 1684 } 1685 } 1686 /* Create and fill new matrix */ 1687 if (scall == MAT_REUSE_MATRIX) { 1688 PetscTruth equal; 1689 1690 c = (Mat_SeqAIJ *)((*B)->data); 1691 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1692 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1693 if (!equal) { 1694 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1695 } 1696 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(PetscInt));CHKERRQ(ierr); 1697 C = *B; 1698 } else { 1699 ierr = MatCreate(A->comm,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE,&C);CHKERRQ(ierr); 1700 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1701 ierr = MatSeqAIJSetPreallocation(C,0,lens);CHKERRQ(ierr); 1702 } 1703 c = (Mat_SeqAIJ *)(C->data); 1704 for (i=0; i<nrows; i++) { 1705 row = irow[i]; 1706 kstart = ai[row]; 1707 kend = kstart + a->ilen[row]; 1708 mat_i = c->i[i]; 1709 mat_j = c->j + mat_i; 1710 mat_a = c->a + mat_i; 1711 mat_ilen = c->ilen + i; 1712 for (k=kstart; k<kend; k++) { 1713 if ((tcol=smap[a->j[k]])) { 1714 *mat_j++ = tcol - 1; 1715 *mat_a++ = a->a[k]; 1716 (*mat_ilen)++; 1717 1718 } 1719 } 1720 } 1721 /* Free work space */ 1722 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1723 ierr = PetscFree(smap);CHKERRQ(ierr); 1724 ierr = PetscFree(lens);CHKERRQ(ierr); 1725 } 1726 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1727 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1728 1729 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1730 *B = C; 1731 PetscFunctionReturn(0); 1732 } 1733 1734 /* 1735 */ 1736 #undef __FUNCT__ 1737 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1738 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1739 { 1740 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1741 PetscErrorCode ierr; 1742 Mat outA; 1743 PetscTruth row_identity,col_identity; 1744 1745 PetscFunctionBegin; 1746 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1747 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1748 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1749 if (!row_identity || !col_identity) { 1750 SETERRQ(PETSC_ERR_ARG_WRONG,"Row and column permutations must be identity for in-place ILU"); 1751 } 1752 1753 outA = inA; 1754 inA->factor = FACTOR_LU; 1755 a->row = row; 1756 a->col = col; 1757 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1758 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1759 1760 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1761 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1762 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1763 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1764 1765 if (!a->solve_work) { /* this matrix may have been factored before */ 1766 ierr = PetscMalloc((inA->m+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1767 } 1768 1769 if (!a->diag) { 1770 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1771 } 1772 ierr = MatLUFactorNumeric_SeqAIJ(inA,info,&outA);CHKERRQ(ierr); 1773 PetscFunctionReturn(0); 1774 } 1775 1776 #include "petscblaslapack.h" 1777 #undef __FUNCT__ 1778 #define __FUNCT__ "MatScale_SeqAIJ" 1779 PetscErrorCode MatScale_SeqAIJ(const PetscScalar *alpha,Mat inA) 1780 { 1781 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1782 PetscBLASInt bnz = (PetscBLASInt)a->nz,one = 1; 1783 PetscErrorCode ierr; 1784 1785 1786 PetscFunctionBegin; 1787 BLASscal_(&bnz,(PetscScalar*)alpha,a->a,&one); 1788 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1789 PetscFunctionReturn(0); 1790 } 1791 1792 #undef __FUNCT__ 1793 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1794 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1795 { 1796 PetscErrorCode ierr; 1797 PetscInt i; 1798 1799 PetscFunctionBegin; 1800 if (scall == MAT_INITIAL_MATRIX) { 1801 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1802 } 1803 1804 for (i=0; i<n; i++) { 1805 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1806 } 1807 PetscFunctionReturn(0); 1808 } 1809 1810 #undef __FUNCT__ 1811 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1812 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1813 { 1814 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1815 PetscErrorCode ierr; 1816 PetscInt row,i,j,k,l,m,n,*idx,*nidx,isz,val; 1817 PetscInt start,end,*ai,*aj; 1818 PetscBT table; 1819 1820 PetscFunctionBegin; 1821 m = A->m; 1822 ai = a->i; 1823 aj = a->j; 1824 1825 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1826 1827 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1828 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1829 1830 for (i=0; i<is_max; i++) { 1831 /* Initialize the two local arrays */ 1832 isz = 0; 1833 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1834 1835 /* Extract the indices, assume there can be duplicate entries */ 1836 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1837 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1838 1839 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1840 for (j=0; j<n ; ++j){ 1841 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1842 } 1843 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1844 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1845 1846 k = 0; 1847 for (j=0; j<ov; j++){ /* for each overlap */ 1848 n = isz; 1849 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1850 row = nidx[k]; 1851 start = ai[row]; 1852 end = ai[row+1]; 1853 for (l = start; l<end ; l++){ 1854 val = aj[l] ; 1855 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1856 } 1857 } 1858 } 1859 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1860 } 1861 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1862 ierr = PetscFree(nidx);CHKERRQ(ierr); 1863 PetscFunctionReturn(0); 1864 } 1865 1866 /* -------------------------------------------------------------- */ 1867 #undef __FUNCT__ 1868 #define __FUNCT__ "MatPermute_SeqAIJ" 1869 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1870 { 1871 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1872 PetscErrorCode ierr; 1873 PetscInt i,nz,m = A->m,n = A->n,*col; 1874 PetscInt *row,*cnew,j,*lens; 1875 IS icolp,irowp; 1876 PetscInt *cwork; 1877 PetscScalar *vwork; 1878 1879 PetscFunctionBegin; 1880 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1881 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1882 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1883 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1884 1885 /* determine lengths of permuted rows */ 1886 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1887 for (i=0; i<m; i++) { 1888 lens[row[i]] = a->i[i+1] - a->i[i]; 1889 } 1890 ierr = MatCreate(A->comm,m,n,m,n,B);CHKERRQ(ierr); 1891 ierr = MatSetType(*B,A->type_name);CHKERRQ(ierr); 1892 ierr = MatSeqAIJSetPreallocation(*B,0,lens);CHKERRQ(ierr); 1893 ierr = PetscFree(lens);CHKERRQ(ierr); 1894 1895 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 1896 for (i=0; i<m; i++) { 1897 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1898 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1899 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1900 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1901 } 1902 ierr = PetscFree(cnew);CHKERRQ(ierr); 1903 (*B)->assembled = PETSC_FALSE; 1904 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1905 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1906 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1907 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1908 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1909 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1910 PetscFunctionReturn(0); 1911 } 1912 1913 #undef __FUNCT__ 1914 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1915 PetscErrorCode MatPrintHelp_SeqAIJ(Mat A) 1916 { 1917 static PetscTruth called = PETSC_FALSE; 1918 MPI_Comm comm = A->comm; 1919 PetscErrorCode ierr; 1920 1921 PetscFunctionBegin; 1922 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1923 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1924 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1925 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1926 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1927 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1928 ierr = (*PetscHelpPrintf)(comm," -mat_no_compressedrow: Do not use compressedrow\n");CHKERRQ(ierr); 1929 PetscFunctionReturn(0); 1930 } 1931 1932 #undef __FUNCT__ 1933 #define __FUNCT__ "MatCopy_SeqAIJ" 1934 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1935 { 1936 PetscErrorCode ierr; 1937 1938 PetscFunctionBegin; 1939 /* If the two matrices have the same copy implementation, use fast copy. */ 1940 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1941 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1942 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1943 1944 if (a->i[A->m] != b->i[B->m]) { 1945 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1946 } 1947 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 1948 } else { 1949 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1950 } 1951 PetscFunctionReturn(0); 1952 } 1953 1954 #undef __FUNCT__ 1955 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1956 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 1957 { 1958 PetscErrorCode ierr; 1959 1960 PetscFunctionBegin; 1961 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1962 PetscFunctionReturn(0); 1963 } 1964 1965 #undef __FUNCT__ 1966 #define __FUNCT__ "MatGetArray_SeqAIJ" 1967 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 1968 { 1969 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1970 PetscFunctionBegin; 1971 *array = a->a; 1972 PetscFunctionReturn(0); 1973 } 1974 1975 #undef __FUNCT__ 1976 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1977 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 1978 { 1979 PetscFunctionBegin; 1980 PetscFunctionReturn(0); 1981 } 1982 1983 #undef __FUNCT__ 1984 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1985 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1986 { 1987 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 1988 PetscErrorCode ierr; 1989 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1990 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1991 PetscScalar *vscale_array; 1992 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1993 Vec w1,w2,w3; 1994 void *fctx = coloring->fctx; 1995 PetscTruth flg; 1996 1997 PetscFunctionBegin; 1998 if (!coloring->w1) { 1999 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2000 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2001 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2002 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2003 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2004 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2005 } 2006 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2007 2008 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2009 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 2010 if (flg) { 2011 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 2012 } else { 2013 PetscTruth assembled; 2014 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2015 if (assembled) { 2016 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2017 } 2018 } 2019 2020 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2021 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2022 2023 /* 2024 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2025 coloring->F for the coarser grids from the finest 2026 */ 2027 if (coloring->F) { 2028 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2029 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2030 if (m1 != m2) { 2031 coloring->F = 0; 2032 } 2033 } 2034 2035 if (coloring->F) { 2036 w1 = coloring->F; 2037 coloring->F = 0; 2038 } else { 2039 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2040 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2041 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2042 } 2043 2044 /* 2045 Compute all the scale factors and share with other processors 2046 */ 2047 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2048 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2049 for (k=0; k<coloring->ncolors; k++) { 2050 /* 2051 Loop over each column associated with color adding the 2052 perturbation to the vector w3. 2053 */ 2054 for (l=0; l<coloring->ncolumns[k]; l++) { 2055 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2056 dx = xx[col]; 2057 if (dx == 0.0) dx = 1.0; 2058 #if !defined(PETSC_USE_COMPLEX) 2059 if (dx < umin && dx >= 0.0) dx = umin; 2060 else if (dx < 0.0 && dx > -umin) dx = -umin; 2061 #else 2062 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2063 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2064 #endif 2065 dx *= epsilon; 2066 vscale_array[col] = 1.0/dx; 2067 } 2068 } 2069 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2070 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2071 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2072 2073 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2074 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2075 2076 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2077 else vscaleforrow = coloring->columnsforrow; 2078 2079 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2080 /* 2081 Loop over each color 2082 */ 2083 for (k=0; k<coloring->ncolors; k++) { 2084 coloring->currentcolor = k; 2085 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2086 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2087 /* 2088 Loop over each column associated with color adding the 2089 perturbation to the vector w3. 2090 */ 2091 for (l=0; l<coloring->ncolumns[k]; l++) { 2092 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2093 dx = xx[col]; 2094 if (dx == 0.0) dx = 1.0; 2095 #if !defined(PETSC_USE_COMPLEX) 2096 if (dx < umin && dx >= 0.0) dx = umin; 2097 else if (dx < 0.0 && dx > -umin) dx = -umin; 2098 #else 2099 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2100 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2101 #endif 2102 dx *= epsilon; 2103 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2104 w3_array[col] += dx; 2105 } 2106 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2107 2108 /* 2109 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2110 */ 2111 2112 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2113 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2114 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2115 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 2116 2117 /* 2118 Loop over rows of vector, putting results into Jacobian matrix 2119 */ 2120 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2121 for (l=0; l<coloring->nrows[k]; l++) { 2122 row = coloring->rows[k][l]; 2123 col = coloring->columnsforrow[k][l]; 2124 y[row] *= vscale_array[vscaleforrow[k][l]]; 2125 srow = row + start; 2126 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2127 } 2128 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2129 } 2130 coloring->currentcolor = k; 2131 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2132 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2133 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2134 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2135 PetscFunctionReturn(0); 2136 } 2137 2138 #include "petscblaslapack.h" 2139 #undef __FUNCT__ 2140 #define __FUNCT__ "MatAXPY_SeqAIJ" 2141 PetscErrorCode MatAXPY_SeqAIJ(const PetscScalar a[],Mat X,Mat Y,MatStructure str) 2142 { 2143 PetscErrorCode ierr; 2144 PetscInt i; 2145 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2146 PetscBLASInt one=1,bnz = (PetscBLASInt)x->nz; 2147 2148 PetscFunctionBegin; 2149 if (str == SAME_NONZERO_PATTERN) { 2150 BLASaxpy_(&bnz,(PetscScalar*)a,x->a,&one,y->a,&one); 2151 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2152 if (y->xtoy && y->XtoY != X) { 2153 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2154 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2155 } 2156 if (!y->xtoy) { /* get xtoy */ 2157 ierr = MatAXPYGetxtoy_Private(X->m,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2158 y->XtoY = X; 2159 } 2160 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += (*a)*(x->a[i]); 2161 PetscLogInfo(0,"MatAXPY_SeqAIJ: ratio of nnz(X)/nnz(Y): %d/%d = %g\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz); 2162 } else { 2163 ierr = MatAXPY_Basic(a,X,Y,str);CHKERRQ(ierr); 2164 } 2165 PetscFunctionReturn(0); 2166 } 2167 2168 #undef __FUNCT__ 2169 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2170 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2171 { 2172 PetscFunctionBegin; 2173 PetscFunctionReturn(0); 2174 } 2175 2176 /* -------------------------------------------------------------------*/ 2177 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2178 MatGetRow_SeqAIJ, 2179 MatRestoreRow_SeqAIJ, 2180 MatMult_SeqAIJ, 2181 /* 4*/ MatMultAdd_SeqAIJ, 2182 MatMultTranspose_SeqAIJ, 2183 MatMultTransposeAdd_SeqAIJ, 2184 MatSolve_SeqAIJ, 2185 MatSolveAdd_SeqAIJ, 2186 MatSolveTranspose_SeqAIJ, 2187 /*10*/ MatSolveTransposeAdd_SeqAIJ, 2188 MatLUFactor_SeqAIJ, 2189 0, 2190 MatRelax_SeqAIJ, 2191 MatTranspose_SeqAIJ, 2192 /*15*/ MatGetInfo_SeqAIJ, 2193 MatEqual_SeqAIJ, 2194 MatGetDiagonal_SeqAIJ, 2195 MatDiagonalScale_SeqAIJ, 2196 MatNorm_SeqAIJ, 2197 /*20*/ 0, 2198 MatAssemblyEnd_SeqAIJ, 2199 MatCompress_SeqAIJ, 2200 MatSetOption_SeqAIJ, 2201 MatZeroEntries_SeqAIJ, 2202 /*25*/ MatZeroRows_SeqAIJ, 2203 MatLUFactorSymbolic_SeqAIJ, 2204 MatLUFactorNumeric_SeqAIJ, 2205 MatCholeskyFactorSymbolic_SeqAIJ, 2206 MatCholeskyFactorNumeric_SeqAIJ, 2207 /*30*/ MatSetUpPreallocation_SeqAIJ, 2208 MatILUFactorSymbolic_SeqAIJ, 2209 MatICCFactorSymbolic_SeqAIJ, 2210 MatGetArray_SeqAIJ, 2211 MatRestoreArray_SeqAIJ, 2212 /*35*/ MatDuplicate_SeqAIJ, 2213 0, 2214 0, 2215 MatILUFactor_SeqAIJ, 2216 0, 2217 /*40*/ MatAXPY_SeqAIJ, 2218 MatGetSubMatrices_SeqAIJ, 2219 MatIncreaseOverlap_SeqAIJ, 2220 MatGetValues_SeqAIJ, 2221 MatCopy_SeqAIJ, 2222 /*45*/ MatPrintHelp_SeqAIJ, 2223 MatScale_SeqAIJ, 2224 0, 2225 0, 2226 MatILUDTFactor_SeqAIJ, 2227 /*50*/ MatSetBlockSize_SeqAIJ, 2228 MatGetRowIJ_SeqAIJ, 2229 MatRestoreRowIJ_SeqAIJ, 2230 MatGetColumnIJ_SeqAIJ, 2231 MatRestoreColumnIJ_SeqAIJ, 2232 /*55*/ MatFDColoringCreate_SeqAIJ, 2233 0, 2234 0, 2235 MatPermute_SeqAIJ, 2236 0, 2237 /*60*/ 0, 2238 MatDestroy_SeqAIJ, 2239 MatView_SeqAIJ, 2240 MatGetPetscMaps_Petsc, 2241 0, 2242 /*65*/ 0, 2243 0, 2244 0, 2245 0, 2246 0, 2247 /*70*/ 0, 2248 0, 2249 MatSetColoring_SeqAIJ, 2250 #if defined(PETSC_HAVE_ADIC) 2251 MatSetValuesAdic_SeqAIJ, 2252 #else 2253 0, 2254 #endif 2255 MatSetValuesAdifor_SeqAIJ, 2256 /*75*/ MatFDColoringApply_SeqAIJ, 2257 0, 2258 0, 2259 0, 2260 0, 2261 /*80*/ 0, 2262 0, 2263 0, 2264 0, 2265 MatLoad_SeqAIJ, 2266 /*85*/ MatIsSymmetric_SeqAIJ, 2267 0, 2268 0, 2269 0, 2270 0, 2271 /*90*/ MatMatMult_SeqAIJ_SeqAIJ, 2272 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2273 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2274 MatPtAP_Basic, 2275 MatPtAPSymbolic_SeqAIJ, 2276 /*95*/ MatPtAPNumeric_SeqAIJ, 2277 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2278 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2279 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2280 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2281 /*100*/MatPtAPNumeric_SeqAIJ_SeqAIJ, 2282 0, 2283 0, 2284 }; 2285 2286 EXTERN_C_BEGIN 2287 #undef __FUNCT__ 2288 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2289 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2290 { 2291 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2292 PetscInt i,nz,n; 2293 2294 PetscFunctionBegin; 2295 2296 nz = aij->maxnz; 2297 n = mat->n; 2298 for (i=0; i<nz; i++) { 2299 aij->j[i] = indices[i]; 2300 } 2301 aij->nz = nz; 2302 for (i=0; i<n; i++) { 2303 aij->ilen[i] = aij->imax[i]; 2304 } 2305 2306 PetscFunctionReturn(0); 2307 } 2308 EXTERN_C_END 2309 2310 #undef __FUNCT__ 2311 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2312 /*@ 2313 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2314 in the matrix. 2315 2316 Input Parameters: 2317 + mat - the SeqAIJ matrix 2318 - indices - the column indices 2319 2320 Level: advanced 2321 2322 Notes: 2323 This can be called if you have precomputed the nonzero structure of the 2324 matrix and want to provide it to the matrix object to improve the performance 2325 of the MatSetValues() operation. 2326 2327 You MUST have set the correct numbers of nonzeros per row in the call to 2328 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2329 2330 MUST be called before any calls to MatSetValues(); 2331 2332 The indices should start with zero, not one. 2333 2334 @*/ 2335 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2336 { 2337 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2338 2339 PetscFunctionBegin; 2340 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2341 PetscValidPointer(indices,2); 2342 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2343 if (f) { 2344 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2345 } else { 2346 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2347 } 2348 PetscFunctionReturn(0); 2349 } 2350 2351 /* ----------------------------------------------------------------------------------------*/ 2352 2353 EXTERN_C_BEGIN 2354 #undef __FUNCT__ 2355 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2356 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 2357 { 2358 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2359 PetscErrorCode ierr; 2360 size_t nz = aij->i[mat->m]; 2361 2362 PetscFunctionBegin; 2363 if (aij->nonew != 1) { 2364 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2365 } 2366 2367 /* allocate space for values if not already there */ 2368 if (!aij->saved_values) { 2369 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2370 } 2371 2372 /* copy values over */ 2373 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2374 PetscFunctionReturn(0); 2375 } 2376 EXTERN_C_END 2377 2378 #undef __FUNCT__ 2379 #define __FUNCT__ "MatStoreValues" 2380 /*@ 2381 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2382 example, reuse of the linear part of a Jacobian, while recomputing the 2383 nonlinear portion. 2384 2385 Collect on Mat 2386 2387 Input Parameters: 2388 . mat - the matrix (currently only AIJ matrices support this option) 2389 2390 Level: advanced 2391 2392 Common Usage, with SNESSolve(): 2393 $ Create Jacobian matrix 2394 $ Set linear terms into matrix 2395 $ Apply boundary conditions to matrix, at this time matrix must have 2396 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2397 $ boundary conditions again will not change the nonzero structure 2398 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2399 $ ierr = MatStoreValues(mat); 2400 $ Call SNESSetJacobian() with matrix 2401 $ In your Jacobian routine 2402 $ ierr = MatRetrieveValues(mat); 2403 $ Set nonlinear terms in matrix 2404 2405 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2406 $ // build linear portion of Jacobian 2407 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2408 $ ierr = MatStoreValues(mat); 2409 $ loop over nonlinear iterations 2410 $ ierr = MatRetrieveValues(mat); 2411 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2412 $ // call MatAssemblyBegin/End() on matrix 2413 $ Solve linear system with Jacobian 2414 $ endloop 2415 2416 Notes: 2417 Matrix must already be assemblied before calling this routine 2418 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2419 calling this routine. 2420 2421 .seealso: MatRetrieveValues() 2422 2423 @*/ 2424 PetscErrorCode MatStoreValues(Mat mat) 2425 { 2426 PetscErrorCode ierr,(*f)(Mat); 2427 2428 PetscFunctionBegin; 2429 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2430 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2431 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2432 2433 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2434 if (f) { 2435 ierr = (*f)(mat);CHKERRQ(ierr); 2436 } else { 2437 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2438 } 2439 PetscFunctionReturn(0); 2440 } 2441 2442 EXTERN_C_BEGIN 2443 #undef __FUNCT__ 2444 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2445 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 2446 { 2447 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2448 PetscErrorCode ierr; 2449 PetscInt nz = aij->i[mat->m]; 2450 2451 PetscFunctionBegin; 2452 if (aij->nonew != 1) { 2453 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2454 } 2455 if (!aij->saved_values) { 2456 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2457 } 2458 /* copy values over */ 2459 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2460 PetscFunctionReturn(0); 2461 } 2462 EXTERN_C_END 2463 2464 #undef __FUNCT__ 2465 #define __FUNCT__ "MatRetrieveValues" 2466 /*@ 2467 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2468 example, reuse of the linear part of a Jacobian, while recomputing the 2469 nonlinear portion. 2470 2471 Collect on Mat 2472 2473 Input Parameters: 2474 . mat - the matrix (currently on AIJ matrices support this option) 2475 2476 Level: advanced 2477 2478 .seealso: MatStoreValues() 2479 2480 @*/ 2481 PetscErrorCode MatRetrieveValues(Mat mat) 2482 { 2483 PetscErrorCode ierr,(*f)(Mat); 2484 2485 PetscFunctionBegin; 2486 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2487 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2488 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2489 2490 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2491 if (f) { 2492 ierr = (*f)(mat);CHKERRQ(ierr); 2493 } else { 2494 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2495 } 2496 PetscFunctionReturn(0); 2497 } 2498 2499 2500 /* --------------------------------------------------------------------------------*/ 2501 #undef __FUNCT__ 2502 #define __FUNCT__ "MatCreateSeqAIJ" 2503 /*@C 2504 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2505 (the default parallel PETSc format). For good matrix assembly performance 2506 the user should preallocate the matrix storage by setting the parameter nz 2507 (or the array nnz). By setting these parameters accurately, performance 2508 during matrix assembly can be increased by more than a factor of 50. 2509 2510 Collective on MPI_Comm 2511 2512 Input Parameters: 2513 + comm - MPI communicator, set to PETSC_COMM_SELF 2514 . m - number of rows 2515 . n - number of columns 2516 . nz - number of nonzeros per row (same for all rows) 2517 - nnz - array containing the number of nonzeros in the various rows 2518 (possibly different for each row) or PETSC_NULL 2519 2520 Output Parameter: 2521 . A - the matrix 2522 2523 Notes: 2524 If nnz is given then nz is ignored 2525 2526 The AIJ format (also called the Yale sparse matrix format or 2527 compressed row storage), is fully compatible with standard Fortran 77 2528 storage. That is, the stored row and column indices can begin at 2529 either one (as in Fortran) or zero. See the users' manual for details. 2530 2531 Specify the preallocated storage with either nz or nnz (not both). 2532 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2533 allocation. For large problems you MUST preallocate memory or you 2534 will get TERRIBLE performance, see the users' manual chapter on matrices. 2535 2536 By default, this format uses inodes (identical nodes) when possible, to 2537 improve numerical efficiency of matrix-vector products and solves. We 2538 search for consecutive rows with the same nonzero structure, thereby 2539 reusing matrix information to achieve increased efficiency. 2540 2541 Options Database Keys: 2542 + -mat_aij_no_inode - Do not use inodes 2543 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2544 - -mat_aij_oneindex - Internally use indexing starting at 1 2545 rather than 0. Note that when calling MatSetValues(), 2546 the user still MUST index entries starting at 0! 2547 2548 Level: intermediate 2549 2550 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2551 2552 @*/ 2553 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2554 { 2555 PetscErrorCode ierr; 2556 2557 PetscFunctionBegin; 2558 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2559 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2560 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2561 PetscFunctionReturn(0); 2562 } 2563 2564 #define SKIP_ALLOCATION -4 2565 2566 #undef __FUNCT__ 2567 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2568 /*@C 2569 MatSeqAIJSetPreallocation - For good matrix assembly performance 2570 the user should preallocate the matrix storage by setting the parameter nz 2571 (or the array nnz). By setting these parameters accurately, performance 2572 during matrix assembly can be increased by more than a factor of 50. 2573 2574 Collective on MPI_Comm 2575 2576 Input Parameters: 2577 + comm - MPI communicator, set to PETSC_COMM_SELF 2578 . m - number of rows 2579 . n - number of columns 2580 . nz - number of nonzeros per row (same for all rows) 2581 - nnz - array containing the number of nonzeros in the various rows 2582 (possibly different for each row) or PETSC_NULL 2583 2584 Output Parameter: 2585 . A - the matrix 2586 2587 Notes: 2588 If nnz is given then nz is ignored 2589 2590 The AIJ format (also called the Yale sparse matrix format or 2591 compressed row storage), is fully compatible with standard Fortran 77 2592 storage. That is, the stored row and column indices can begin at 2593 either one (as in Fortran) or zero. See the users' manual for details. 2594 2595 Specify the preallocated storage with either nz or nnz (not both). 2596 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2597 allocation. For large problems you MUST preallocate memory or you 2598 will get TERRIBLE performance, see the users' manual chapter on matrices. 2599 2600 By default, this format uses inodes (identical nodes) when possible, to 2601 improve numerical efficiency of matrix-vector products and solves. We 2602 search for consecutive rows with the same nonzero structure, thereby 2603 reusing matrix information to achieve increased efficiency. 2604 2605 Options Database Keys: 2606 + -mat_aij_no_inode - Do not use inodes 2607 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2608 - -mat_aij_oneindex - Internally use indexing starting at 1 2609 rather than 0. Note that when calling MatSetValues(), 2610 the user still MUST index entries starting at 0! 2611 2612 Level: intermediate 2613 2614 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2615 2616 @*/ 2617 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2618 { 2619 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2620 2621 PetscFunctionBegin; 2622 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2623 if (f) { 2624 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2625 } 2626 PetscFunctionReturn(0); 2627 } 2628 2629 EXTERN_C_BEGIN 2630 #undef __FUNCT__ 2631 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2632 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2633 { 2634 Mat_SeqAIJ *b; 2635 size_t len = 0; 2636 PetscTruth skipallocation = PETSC_FALSE; 2637 PetscErrorCode ierr; 2638 PetscInt i; 2639 2640 PetscFunctionBegin; 2641 2642 if (nz == SKIP_ALLOCATION) { 2643 skipallocation = PETSC_TRUE; 2644 nz = 0; 2645 } 2646 2647 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2648 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2649 if (nnz) { 2650 for (i=0; i<B->m; i++) { 2651 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2652 if (nnz[i] > B->n) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->n); 2653 } 2654 } 2655 2656 B->preallocated = PETSC_TRUE; 2657 b = (Mat_SeqAIJ*)B->data; 2658 2659 ierr = PetscMalloc((B->m+1)*sizeof(PetscInt),&b->imax);CHKERRQ(ierr); 2660 if (!nnz) { 2661 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2662 else if (nz <= 0) nz = 1; 2663 for (i=0; i<B->m; i++) b->imax[i] = nz; 2664 nz = nz*B->m; 2665 } else { 2666 nz = 0; 2667 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2668 } 2669 2670 if (!skipallocation) { 2671 /* allocate the matrix space */ 2672 len = ((size_t) nz)*(sizeof(PetscInt) + sizeof(PetscScalar)) + (B->m+1)*sizeof(PetscInt); 2673 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2674 b->j = (PetscInt*)(b->a + nz); 2675 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2676 b->i = b->j + nz; 2677 b->i[0] = 0; 2678 for (i=1; i<B->m+1; i++) { 2679 b->i[i] = b->i[i-1] + b->imax[i-1]; 2680 } 2681 b->singlemalloc = PETSC_TRUE; 2682 b->freedata = PETSC_TRUE; 2683 } else { 2684 b->freedata = PETSC_FALSE; 2685 } 2686 2687 /* b->ilen will count nonzeros in each row so far. */ 2688 ierr = PetscMalloc((B->m+1)*sizeof(PetscInt),&b->ilen);CHKERRQ(ierr); 2689 ierr = PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(PetscInt)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2690 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2691 2692 b->nz = 0; 2693 b->maxnz = nz; 2694 B->info.nz_unneeded = (double)b->maxnz; 2695 PetscFunctionReturn(0); 2696 } 2697 EXTERN_C_END 2698 2699 /*MC 2700 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 2701 based on compressed sparse row format. 2702 2703 Options Database Keys: 2704 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 2705 2706 Level: beginner 2707 2708 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 2709 M*/ 2710 2711 EXTERN_C_BEGIN 2712 #undef __FUNCT__ 2713 #define __FUNCT__ "MatCreate_SeqAIJ" 2714 PetscErrorCode MatCreate_SeqAIJ(Mat B) 2715 { 2716 Mat_SeqAIJ *b; 2717 PetscErrorCode ierr; 2718 PetscMPIInt size; 2719 2720 PetscFunctionBegin; 2721 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2722 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2723 2724 B->m = B->M = PetscMax(B->m,B->M); 2725 B->n = B->N = PetscMax(B->n,B->N); 2726 2727 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2728 B->data = (void*)b; 2729 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2730 B->factor = 0; 2731 B->lupivotthreshold = 1.0; 2732 B->mapping = 0; 2733 ierr = PetscOptionsGetReal(B->prefix,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2734 ierr = PetscOptionsHasName(B->prefix,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2735 b->row = 0; 2736 b->col = 0; 2737 b->icol = 0; 2738 b->reallocs = 0; 2739 b->lu_shift = PETSC_FALSE; 2740 2741 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2742 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2743 2744 b->sorted = PETSC_FALSE; 2745 b->ignorezeroentries = PETSC_FALSE; 2746 b->roworiented = PETSC_TRUE; 2747 b->nonew = 0; 2748 b->diag = 0; 2749 b->solve_work = 0; 2750 B->spptr = 0; 2751 b->inode.use = PETSC_TRUE; 2752 b->inode.node_count = 0; 2753 b->inode.size = 0; 2754 b->inode.limit = 5; 2755 b->inode.max_limit = 5; 2756 b->saved_values = 0; 2757 b->idiag = 0; 2758 b->ssor = 0; 2759 b->keepzeroedrows = PETSC_FALSE; 2760 b->xtoy = 0; 2761 b->XtoY = 0; 2762 b->compressedrow.use = PETSC_FALSE; 2763 b->compressedrow.nrows = B->m; 2764 b->compressedrow.i = PETSC_NULL; 2765 b->compressedrow.rindex = PETSC_NULL; 2766 b->compressedrow.checked = PETSC_FALSE; 2767 B->same_nonzero = PETSC_FALSE; 2768 2769 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2770 2771 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2772 "MatSeqAIJSetColumnIndices_SeqAIJ", 2773 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2774 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2775 "MatStoreValues_SeqAIJ", 2776 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2777 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2778 "MatRetrieveValues_SeqAIJ", 2779 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2780 #if !defined(PETSC_USE_64BIT_INT) 2781 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2782 "MatConvert_SeqAIJ_SeqSBAIJ", 2783 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2784 #endif 2785 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2786 "MatConvert_SeqAIJ_SeqBAIJ", 2787 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2788 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 2789 "MatIsTranspose_SeqAIJ", 2790 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 2791 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2792 "MatSeqAIJSetPreallocation_SeqAIJ", 2793 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2794 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 2795 "MatReorderForNonzeroDiagonal_SeqAIJ", 2796 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 2797 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatAdjustForInodes_C", 2798 "MatAdjustForInodes_SeqAIJ", 2799 MatAdjustForInodes_SeqAIJ);CHKERRQ(ierr); 2800 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJGetInodeSizes_C", 2801 "MatSeqAIJGetInodeSizes_SeqAIJ", 2802 MatSeqAIJGetInodeSizes_SeqAIJ);CHKERRQ(ierr); 2803 PetscFunctionReturn(0); 2804 } 2805 EXTERN_C_END 2806 2807 #undef __FUNCT__ 2808 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2809 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2810 { 2811 Mat C; 2812 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2813 PetscErrorCode ierr; 2814 PetscInt i,m = A->m; 2815 size_t len; 2816 2817 PetscFunctionBegin; 2818 *B = 0; 2819 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2820 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 2821 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2822 2823 c = (Mat_SeqAIJ*)C->data; 2824 2825 C->factor = A->factor; 2826 C->lupivotthreshold = A->lupivotthreshold; 2827 2828 c->row = 0; 2829 c->col = 0; 2830 c->icol = 0; 2831 c->reallocs = 0; 2832 c->lu_shift = PETSC_FALSE; 2833 2834 C->assembled = PETSC_TRUE; 2835 2836 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->imax);CHKERRQ(ierr); 2837 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->ilen);CHKERRQ(ierr); 2838 for (i=0; i<m; i++) { 2839 c->imax[i] = a->imax[i]; 2840 c->ilen[i] = a->ilen[i]; 2841 } 2842 2843 /* allocate the matrix space */ 2844 c->singlemalloc = PETSC_TRUE; 2845 len = ((size_t) (m+1))*sizeof(PetscInt)+(a->i[m])*(sizeof(PetscScalar)+sizeof(PetscInt)); 2846 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2847 c->j = (PetscInt*)(c->a + a->i[m] ); 2848 c->i = c->j + a->i[m]; 2849 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2850 if (m > 0) { 2851 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 2852 if (cpvalues == MAT_COPY_VALUES) { 2853 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2854 } else { 2855 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2856 } 2857 } 2858 2859 ierr = PetscLogObjectMemory(C,len+2*(m+1)*sizeof(PetscInt)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2860 c->sorted = a->sorted; 2861 c->ignorezeroentries = a->ignorezeroentries; 2862 c->roworiented = a->roworiented; 2863 c->nonew = a->nonew; 2864 if (a->diag) { 2865 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2866 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2867 for (i=0; i<m; i++) { 2868 c->diag[i] = a->diag[i]; 2869 } 2870 } else c->diag = 0; 2871 c->solve_work = 0; 2872 c->inode.use = a->inode.use; 2873 c->inode.limit = a->inode.limit; 2874 c->inode.max_limit = a->inode.max_limit; 2875 if (a->inode.size){ 2876 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->inode.size);CHKERRQ(ierr); 2877 c->inode.node_count = a->inode.node_count; 2878 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2879 } else { 2880 c->inode.size = 0; 2881 c->inode.node_count = 0; 2882 } 2883 c->saved_values = 0; 2884 c->idiag = 0; 2885 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2886 c->ssor = 0; 2887 c->keepzeroedrows = a->keepzeroedrows; 2888 c->freedata = PETSC_TRUE; 2889 c->xtoy = 0; 2890 c->XtoY = 0; 2891 2892 c->nz = a->nz; 2893 c->maxnz = a->maxnz; 2894 C->preallocated = PETSC_TRUE; 2895 2896 c->compressedrow.use = a->compressedrow.use; 2897 c->compressedrow.nrows = a->compressedrow.nrows; 2898 c->compressedrow.checked = a->compressedrow.checked; 2899 if ( a->compressedrow.checked && a->compressedrow.use){ 2900 i = a->compressedrow.nrows; 2901 ierr = PetscMalloc((2*i+1)*sizeof(PetscInt),&c->compressedrow.i);CHKERRQ(ierr); 2902 c->compressedrow.rindex = c->compressedrow.i + i + 1; 2903 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 2904 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 2905 } else { 2906 c->compressedrow.use = PETSC_FALSE; 2907 c->compressedrow.i = PETSC_NULL; 2908 c->compressedrow.rindex = PETSC_NULL; 2909 } 2910 C->same_nonzero = A->same_nonzero; 2911 2912 *B = C; 2913 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2914 PetscFunctionReturn(0); 2915 } 2916 2917 #undef __FUNCT__ 2918 #define __FUNCT__ "MatLoad_SeqAIJ" 2919 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer,const MatType type,Mat *A) 2920 { 2921 Mat_SeqAIJ *a; 2922 Mat B; 2923 PetscErrorCode ierr; 2924 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 2925 int fd; 2926 PetscMPIInt size; 2927 MPI_Comm comm; 2928 2929 PetscFunctionBegin; 2930 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2931 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2932 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2933 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2934 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2935 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2936 M = header[1]; N = header[2]; nz = header[3]; 2937 2938 if (nz < 0) { 2939 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2940 } 2941 2942 /* read in row lengths */ 2943 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2944 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2945 2946 /* check if sum of rowlengths is same as nz */ 2947 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 2948 if (sum != nz) SETERRQ2(PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 2949 2950 /* create our matrix */ 2951 ierr = MatCreate(comm,PETSC_DECIDE,PETSC_DECIDE,M,N,&B);CHKERRQ(ierr); 2952 ierr = MatSetType(B,type);CHKERRQ(ierr); 2953 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 2954 a = (Mat_SeqAIJ*)B->data; 2955 2956 /* read in column indices and adjust for Fortran indexing*/ 2957 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2958 2959 /* read in nonzero values */ 2960 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2961 2962 /* set matrix "i" values */ 2963 a->i[0] = 0; 2964 for (i=1; i<= M; i++) { 2965 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2966 a->ilen[i-1] = rowlengths[i-1]; 2967 } 2968 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2969 2970 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2971 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2972 *A = B; 2973 PetscFunctionReturn(0); 2974 } 2975 2976 #undef __FUNCT__ 2977 #define __FUNCT__ "MatEqual_SeqAIJ" 2978 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2979 { 2980 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2981 PetscErrorCode ierr; 2982 2983 PetscFunctionBegin; 2984 /* If the matrix dimensions are not equal,or no of nonzeros */ 2985 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)) { 2986 *flg = PETSC_FALSE; 2987 PetscFunctionReturn(0); 2988 } 2989 2990 /* if the a->i are the same */ 2991 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 2992 if (!*flg) PetscFunctionReturn(0); 2993 2994 /* if a->j are the same */ 2995 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 2996 if (!*flg) PetscFunctionReturn(0); 2997 2998 /* if a->a are the same */ 2999 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3000 3001 PetscFunctionReturn(0); 3002 3003 } 3004 3005 #undef __FUNCT__ 3006 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3007 /*@C 3008 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3009 provided by the user. 3010 3011 Coolective on MPI_Comm 3012 3013 Input Parameters: 3014 + comm - must be an MPI communicator of size 1 3015 . m - number of rows 3016 . n - number of columns 3017 . i - row indices 3018 . j - column indices 3019 - a - matrix values 3020 3021 Output Parameter: 3022 . mat - the matrix 3023 3024 Level: intermediate 3025 3026 Notes: 3027 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3028 once the matrix is destroyed 3029 3030 You cannot set new nonzero locations into this matrix, that will generate an error. 3031 3032 The i and j indices are 0 based 3033 3034 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 3035 3036 @*/ 3037 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3038 { 3039 PetscErrorCode ierr; 3040 PetscInt ii; 3041 Mat_SeqAIJ *aij; 3042 3043 PetscFunctionBegin; 3044 ierr = MatCreate(comm,m,n,m,n,mat);CHKERRQ(ierr); 3045 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3046 ierr = MatSeqAIJSetPreallocation(*mat,SKIP_ALLOCATION,0);CHKERRQ(ierr); 3047 aij = (Mat_SeqAIJ*)(*mat)->data; 3048 3049 if (i[0] != 0) { 3050 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3051 } 3052 aij->i = i; 3053 aij->j = j; 3054 aij->a = a; 3055 aij->singlemalloc = PETSC_FALSE; 3056 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3057 aij->freedata = PETSC_FALSE; 3058 3059 for (ii=0; ii<m; ii++) { 3060 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3061 #if defined(PETSC_USE_DEBUG) 3062 if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 3063 #endif 3064 } 3065 #if defined(PETSC_USE_DEBUG) 3066 for (ii=0; ii<aij->i[m]; ii++) { 3067 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3068 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3069 } 3070 #endif 3071 3072 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3073 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3074 PetscFunctionReturn(0); 3075 } 3076 3077 #undef __FUNCT__ 3078 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3079 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3080 { 3081 PetscErrorCode ierr; 3082 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3083 3084 PetscFunctionBegin; 3085 if (coloring->ctype == IS_COLORING_LOCAL) { 3086 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3087 a->coloring = coloring; 3088 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3089 PetscInt i,*larray; 3090 ISColoring ocoloring; 3091 ISColoringValue *colors; 3092 3093 /* set coloring for diagonal portion */ 3094 ierr = PetscMalloc((A->n+1)*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3095 for (i=0; i<A->n; i++) { 3096 larray[i] = i; 3097 } 3098 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3099 ierr = PetscMalloc((A->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3100 for (i=0; i<A->n; i++) { 3101 colors[i] = coloring->colors[larray[i]]; 3102 } 3103 ierr = PetscFree(larray);CHKERRQ(ierr); 3104 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 3105 a->coloring = ocoloring; 3106 } 3107 PetscFunctionReturn(0); 3108 } 3109 3110 #if defined(PETSC_HAVE_ADIC) 3111 EXTERN_C_BEGIN 3112 #include "adic/ad_utils.h" 3113 EXTERN_C_END 3114 3115 #undef __FUNCT__ 3116 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3117 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3118 { 3119 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3120 PetscInt m = A->m,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3121 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3122 ISColoringValue *color; 3123 3124 PetscFunctionBegin; 3125 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3126 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3127 color = a->coloring->colors; 3128 /* loop over rows */ 3129 for (i=0; i<m; i++) { 3130 nz = ii[i+1] - ii[i]; 3131 /* loop over columns putting computed value into matrix */ 3132 for (j=0; j<nz; j++) { 3133 *v++ = values[color[*jj++]]; 3134 } 3135 values += nlen; /* jump to next row of derivatives */ 3136 } 3137 PetscFunctionReturn(0); 3138 } 3139 #endif 3140 3141 #undef __FUNCT__ 3142 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3143 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3144 { 3145 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3146 PetscInt m = A->m,*ii = a->i,*jj = a->j,nz,i,j; 3147 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3148 ISColoringValue *color; 3149 3150 PetscFunctionBegin; 3151 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3152 color = a->coloring->colors; 3153 /* loop over rows */ 3154 for (i=0; i<m; i++) { 3155 nz = ii[i+1] - ii[i]; 3156 /* loop over columns putting computed value into matrix */ 3157 for (j=0; j<nz; j++) { 3158 *v++ = values[color[*jj++]]; 3159 } 3160 values += nl; /* jump to next row of derivatives */ 3161 } 3162 PetscFunctionReturn(0); 3163 } 3164 3165