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