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