1 #define PETSCMAT_DLL 2 3 4 /* 5 Defines the basic matrix operations for the AIJ (compressed row) 6 matrix storage format. 7 */ 8 9 10 #include "../src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 11 #include "petscblaslapack.h" 12 #include "petscbt.h" 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "MatDiagonalSet_SeqAIJ" 16 PetscErrorCode PETSCMAT_DLLEXPORT MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 17 { 18 PetscErrorCode ierr; 19 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 20 PetscInt i,*diag, m = Y->rmap->n; 21 MatScalar *aa = aij->a; 22 PetscScalar *v; 23 PetscTruth missing; 24 25 PetscFunctionBegin; 26 if (Y->assembled) { 27 ierr = MatMissingDiagonal_SeqAIJ(Y,&missing,PETSC_NULL);CHKERRQ(ierr); 28 if (!missing) { 29 diag = aij->diag; 30 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 31 if (is == INSERT_VALUES) { 32 for (i=0; i<m; i++) { 33 aa[diag[i]] = v[i]; 34 } 35 } else { 36 for (i=0; i<m; i++) { 37 aa[diag[i]] += v[i]; 38 } 39 } 40 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 41 PetscFunctionReturn(0); 42 } 43 } 44 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 45 PetscFunctionReturn(0); 46 } 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 50 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *m,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 51 { 52 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 53 PetscErrorCode ierr; 54 PetscInt i,ishift; 55 56 PetscFunctionBegin; 57 *m = A->rmap->n; 58 if (!ia) PetscFunctionReturn(0); 59 ishift = 0; 60 if (symmetric && !A->structurally_symmetric) { 61 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 62 } else if (oshift == 1) { 63 PetscInt nz = a->i[A->rmap->n]; 64 /* malloc space and add 1 to i and j indices */ 65 ierr = PetscMalloc((A->rmap->n+1)*sizeof(PetscInt),ia);CHKERRQ(ierr); 66 for (i=0; i<A->rmap->n+1; i++) (*ia)[i] = a->i[i] + 1; 67 if (ja) { 68 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),ja);CHKERRQ(ierr); 69 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 70 } 71 } else { 72 *ia = a->i; 73 if (ja) *ja = a->j; 74 } 75 PetscFunctionReturn(0); 76 } 77 78 #undef __FUNCT__ 79 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 80 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 81 { 82 PetscErrorCode ierr; 83 84 PetscFunctionBegin; 85 if (!ia) PetscFunctionReturn(0); 86 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 87 ierr = PetscFree(*ia);CHKERRQ(ierr); 88 if (ja) {ierr = PetscFree(*ja);CHKERRQ(ierr);} 89 } 90 PetscFunctionReturn(0); 91 } 92 93 #undef __FUNCT__ 94 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 95 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 96 { 97 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 98 PetscErrorCode ierr; 99 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 100 PetscInt nz = a->i[m],row,*jj,mr,col; 101 102 PetscFunctionBegin; 103 *nn = n; 104 if (!ia) PetscFunctionReturn(0); 105 if (symmetric) { 106 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 107 } else { 108 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 109 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 110 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 111 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 112 jj = a->j; 113 for (i=0; i<nz; i++) { 114 collengths[jj[i]]++; 115 } 116 cia[0] = oshift; 117 for (i=0; i<n; i++) { 118 cia[i+1] = cia[i] + collengths[i]; 119 } 120 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 121 jj = a->j; 122 for (row=0; row<m; row++) { 123 mr = a->i[row+1] - a->i[row]; 124 for (i=0; i<mr; i++) { 125 col = *jj++; 126 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 127 } 128 } 129 ierr = PetscFree(collengths);CHKERRQ(ierr); 130 *ia = cia; *ja = cja; 131 } 132 PetscFunctionReturn(0); 133 } 134 135 #undef __FUNCT__ 136 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 137 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 138 { 139 PetscErrorCode ierr; 140 141 PetscFunctionBegin; 142 if (!ia) PetscFunctionReturn(0); 143 144 ierr = PetscFree(*ia);CHKERRQ(ierr); 145 ierr = PetscFree(*ja);CHKERRQ(ierr); 146 147 PetscFunctionReturn(0); 148 } 149 150 #undef __FUNCT__ 151 #define __FUNCT__ "MatSetValuesRow_SeqAIJ" 152 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 153 { 154 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 155 PetscInt *ai = a->i; 156 PetscErrorCode ierr; 157 158 PetscFunctionBegin; 159 ierr = PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "MatSetValues_SeqAIJ" 165 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 166 { 167 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 168 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 169 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 170 PetscErrorCode ierr; 171 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 172 MatScalar *ap,value,*aa = a->a; 173 PetscTruth ignorezeroentries = a->ignorezeroentries; 174 PetscTruth roworiented = a->roworiented; 175 176 PetscFunctionBegin; 177 if (v) PetscValidScalarPointer(v,6); 178 for (k=0; k<m; k++) { /* loop over added rows */ 179 row = im[k]; 180 if (row < 0) continue; 181 #if defined(PETSC_USE_DEBUG) 182 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 183 #endif 184 rp = aj + ai[row]; ap = aa + ai[row]; 185 rmax = imax[row]; nrow = ailen[row]; 186 low = 0; 187 high = nrow; 188 for (l=0; l<n; l++) { /* loop over added columns */ 189 if (in[l] < 0) continue; 190 #if defined(PETSC_USE_DEBUG) 191 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 192 #endif 193 col = in[l]; 194 if (v) { 195 if (roworiented) { 196 value = v[l + k*n]; 197 } else { 198 value = v[k + l*m]; 199 } 200 } else { 201 value = 0.; 202 } 203 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 204 205 if (col <= lastcol) low = 0; else high = nrow; 206 lastcol = col; 207 while (high-low > 5) { 208 t = (low+high)/2; 209 if (rp[t] > col) high = t; 210 else low = t; 211 } 212 for (i=low; i<high; i++) { 213 if (rp[i] > col) break; 214 if (rp[i] == col) { 215 if (is == ADD_VALUES) ap[i] += value; 216 else ap[i] = value; 217 low = i + 1; 218 goto noinsert; 219 } 220 } 221 if (value == 0.0 && ignorezeroentries) goto noinsert; 222 if (nonew == 1) goto noinsert; 223 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 224 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 225 N = nrow++ - 1; a->nz++; high++; 226 /* shift up all the later entries in this row */ 227 for (ii=N; ii>=i; ii--) { 228 rp[ii+1] = rp[ii]; 229 ap[ii+1] = ap[ii]; 230 } 231 rp[i] = col; 232 ap[i] = value; 233 low = i + 1; 234 noinsert:; 235 } 236 ailen[row] = nrow; 237 } 238 A->same_nonzero = PETSC_FALSE; 239 PetscFunctionReturn(0); 240 } 241 242 243 #undef __FUNCT__ 244 #define __FUNCT__ "MatGetValues_SeqAIJ" 245 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 246 { 247 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 248 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 249 PetscInt *ai = a->i,*ailen = a->ilen; 250 MatScalar *ap,*aa = a->a; 251 252 PetscFunctionBegin; 253 for (k=0; k<m; k++) { /* loop over rows */ 254 row = im[k]; 255 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 256 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 257 rp = aj + ai[row]; ap = aa + ai[row]; 258 nrow = ailen[row]; 259 for (l=0; l<n; l++) { /* loop over columns */ 260 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 261 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 262 col = in[l] ; 263 high = nrow; low = 0; /* assume unsorted */ 264 while (high-low > 5) { 265 t = (low+high)/2; 266 if (rp[t] > col) high = t; 267 else low = t; 268 } 269 for (i=low; i<high; i++) { 270 if (rp[i] > col) break; 271 if (rp[i] == col) { 272 *v++ = ap[i]; 273 goto finished; 274 } 275 } 276 *v++ = 0.0; 277 finished:; 278 } 279 } 280 PetscFunctionReturn(0); 281 } 282 283 284 #undef __FUNCT__ 285 #define __FUNCT__ "MatView_SeqAIJ_Binary" 286 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 287 { 288 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 289 PetscErrorCode ierr; 290 PetscInt i,*col_lens; 291 int fd; 292 293 PetscFunctionBegin; 294 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 295 ierr = PetscMalloc((4+A->rmap->n)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 296 col_lens[0] = MAT_FILE_CLASSID; 297 col_lens[1] = A->rmap->n; 298 col_lens[2] = A->cmap->n; 299 col_lens[3] = a->nz; 300 301 /* store lengths of each row and write (including header) to file */ 302 for (i=0; i<A->rmap->n; i++) { 303 col_lens[4+i] = a->i[i+1] - a->i[i]; 304 } 305 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 306 ierr = PetscFree(col_lens);CHKERRQ(ierr); 307 308 /* store column indices (zero start index) */ 309 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 310 311 /* store nonzero values */ 312 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 EXTERN PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 317 318 #undef __FUNCT__ 319 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 320 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 321 { 322 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 323 PetscErrorCode ierr; 324 PetscInt i,j,m = A->rmap->n,shift=0; 325 const char *name; 326 PetscViewerFormat format; 327 328 PetscFunctionBegin; 329 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 330 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 331 if (format == PETSC_VIEWER_ASCII_MATLAB) { 332 PetscInt nofinalvalue = 0; 333 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-!shift)) { 334 nofinalvalue = 1; 335 } 336 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 337 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap->n);CHKERRQ(ierr); 338 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 339 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 340 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 341 342 for (i=0; i<m; i++) { 343 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 344 #if defined(PETSC_USE_COMPLEX) 345 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); 346 #else 347 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 348 #endif 349 } 350 } 351 if (nofinalvalue) { 352 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap->n,0.0);CHKERRQ(ierr); 353 } 354 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 355 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 356 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 357 PetscFunctionReturn(0); 358 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 359 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 360 for (i=0; i<m; i++) { 361 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 362 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 363 #if defined(PETSC_USE_COMPLEX) 364 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 365 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 366 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 367 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 368 } else if (PetscRealPart(a->a[j]) != 0.0) { 369 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 370 } 371 #else 372 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 373 #endif 374 } 375 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 376 } 377 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 378 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 379 PetscInt nzd=0,fshift=1,*sptr; 380 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 381 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&sptr);CHKERRQ(ierr); 382 for (i=0; i<m; i++) { 383 sptr[i] = nzd+1; 384 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 385 if (a->j[j] >= i) { 386 #if defined(PETSC_USE_COMPLEX) 387 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 388 #else 389 if (a->a[j] != 0.0) nzd++; 390 #endif 391 } 392 } 393 } 394 sptr[m] = nzd+1; 395 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 396 for (i=0; i<m+1; i+=6) { 397 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);} 398 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);} 399 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);} 400 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 401 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 402 else {ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr);} 403 } 404 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 405 ierr = PetscFree(sptr);CHKERRQ(ierr); 406 for (i=0; i<m; i++) { 407 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 408 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 409 } 410 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 411 } 412 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 413 for (i=0; i<m; i++) { 414 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 415 if (a->j[j] >= i) { 416 #if defined(PETSC_USE_COMPLEX) 417 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 418 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 419 } 420 #else 421 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 422 #endif 423 } 424 } 425 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 426 } 427 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 428 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 429 PetscInt cnt = 0,jcnt; 430 PetscScalar value; 431 432 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 433 for (i=0; i<m; i++) { 434 jcnt = 0; 435 for (j=0; j<A->cmap->n; j++) { 436 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 437 value = a->a[cnt++]; 438 jcnt++; 439 } else { 440 value = 0.0; 441 } 442 #if defined(PETSC_USE_COMPLEX) 443 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 444 #else 445 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 446 #endif 447 } 448 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 449 } 450 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 451 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 452 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 453 #if defined(PETSC_USE_COMPLEX) 454 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix complex general\n");CHKERRQ(ierr); 455 #else 456 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix real general\n");CHKERRQ(ierr); 457 #endif 458 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %D\n", m, A->cmap->n, a->nz);CHKERRQ(ierr); 459 for (i=0; i<m; i++) { 460 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 461 #if defined(PETSC_USE_COMPLEX) 462 if (PetscImaginaryPart(a->a[j]) > 0.0) { 463 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G %G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 464 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 465 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G -%G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 466 } else { 467 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 468 } 469 #else 470 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %G\n", i+shift, a->j[j]+shift, a->a[j]);CHKERRQ(ierr); 471 #endif 472 } 473 } 474 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 475 } else { 476 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 477 if (A->factortype){ 478 for (i=0; i<m; i++) { 479 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 480 /* L part */ 481 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 482 #if defined(PETSC_USE_COMPLEX) 483 if (PetscImaginaryPart(a->a[j]) > 0.0) { 484 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 485 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 486 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 487 } else { 488 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 489 } 490 #else 491 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 492 #endif 493 } 494 /* diagonal */ 495 j = a->diag[i]; 496 #if defined(PETSC_USE_COMPLEX) 497 if (PetscImaginaryPart(a->a[j]) > 0.0) { 498 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 499 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 500 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 501 } else { 502 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 503 } 504 #else 505 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 506 #endif 507 508 /* U part */ 509 for (j=a->diag[i+1]+1+shift; j<a->diag[i]+shift; j++) { 510 #if defined(PETSC_USE_COMPLEX) 511 if (PetscImaginaryPart(a->a[j]) > 0.0) { 512 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 513 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 514 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 515 } else { 516 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 517 } 518 #else 519 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 520 #endif 521 } 522 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 523 } 524 } else { 525 for (i=0; i<m; i++) { 526 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 527 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 528 #if defined(PETSC_USE_COMPLEX) 529 if (PetscImaginaryPart(a->a[j]) > 0.0) { 530 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 531 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 532 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 533 } else { 534 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 535 } 536 #else 537 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 538 #endif 539 } 540 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 541 } 542 } 543 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 544 } 545 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 546 PetscFunctionReturn(0); 547 } 548 549 #undef __FUNCT__ 550 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 551 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 552 { 553 Mat A = (Mat) Aa; 554 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 555 PetscErrorCode ierr; 556 PetscInt i,j,m = A->rmap->n,color; 557 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 558 PetscViewer viewer; 559 PetscViewerFormat format; 560 561 PetscFunctionBegin; 562 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 563 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 564 565 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 566 /* loop over matrix elements drawing boxes */ 567 568 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 569 /* Blue for negative, Cyan for zero and Red for positive */ 570 color = PETSC_DRAW_BLUE; 571 for (i=0; i<m; i++) { 572 y_l = m - i - 1.0; y_r = y_l + 1.0; 573 for (j=a->i[i]; j<a->i[i+1]; j++) { 574 x_l = a->j[j] ; x_r = x_l + 1.0; 575 #if defined(PETSC_USE_COMPLEX) 576 if (PetscRealPart(a->a[j]) >= 0.) continue; 577 #else 578 if (a->a[j] >= 0.) continue; 579 #endif 580 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 581 } 582 } 583 color = PETSC_DRAW_CYAN; 584 for (i=0; i<m; i++) { 585 y_l = m - i - 1.0; y_r = y_l + 1.0; 586 for (j=a->i[i]; j<a->i[i+1]; j++) { 587 x_l = a->j[j]; x_r = x_l + 1.0; 588 if (a->a[j] != 0.) continue; 589 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 590 } 591 } 592 color = PETSC_DRAW_RED; 593 for (i=0; i<m; i++) { 594 y_l = m - i - 1.0; y_r = y_l + 1.0; 595 for (j=a->i[i]; j<a->i[i+1]; j++) { 596 x_l = a->j[j]; x_r = x_l + 1.0; 597 #if defined(PETSC_USE_COMPLEX) 598 if (PetscRealPart(a->a[j]) <= 0.) continue; 599 #else 600 if (a->a[j] <= 0.) continue; 601 #endif 602 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 603 } 604 } 605 } else { 606 /* use contour shading to indicate magnitude of values */ 607 /* first determine max of all nonzero values */ 608 PetscInt nz = a->nz,count; 609 PetscDraw popup; 610 PetscReal scale; 611 612 for (i=0; i<nz; i++) { 613 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 614 } 615 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 616 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 617 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 618 count = 0; 619 for (i=0; i<m; i++) { 620 y_l = m - i - 1.0; y_r = y_l + 1.0; 621 for (j=a->i[i]; j<a->i[i+1]; j++) { 622 x_l = a->j[j]; x_r = x_l + 1.0; 623 color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count])); 624 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 625 count++; 626 } 627 } 628 } 629 PetscFunctionReturn(0); 630 } 631 632 #undef __FUNCT__ 633 #define __FUNCT__ "MatView_SeqAIJ_Draw" 634 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 635 { 636 PetscErrorCode ierr; 637 PetscDraw draw; 638 PetscReal xr,yr,xl,yl,h,w; 639 PetscTruth isnull; 640 641 PetscFunctionBegin; 642 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 643 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 644 if (isnull) PetscFunctionReturn(0); 645 646 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 647 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 648 xr += w; yr += h; xl = -w; yl = -h; 649 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 650 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 651 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 652 PetscFunctionReturn(0); 653 } 654 655 #undef __FUNCT__ 656 #define __FUNCT__ "MatView_SeqAIJ" 657 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 658 { 659 PetscErrorCode ierr; 660 PetscTruth iascii,isbinary,isdraw; 661 662 PetscFunctionBegin; 663 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 664 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 665 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 666 if (iascii) { 667 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 668 } else if (isbinary) { 669 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 670 } else if (isdraw) { 671 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 672 } else { 673 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 674 } 675 ierr = MatView_SeqAIJ_Inode(A,viewer);CHKERRQ(ierr); 676 PetscFunctionReturn(0); 677 } 678 679 #undef __FUNCT__ 680 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 681 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 682 { 683 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 684 PetscErrorCode ierr; 685 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 686 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 687 MatScalar *aa = a->a,*ap; 688 PetscReal ratio=0.6; 689 690 PetscFunctionBegin; 691 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 692 693 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 694 for (i=1; i<m; i++) { 695 /* move each row back by the amount of empty slots (fshift) before it*/ 696 fshift += imax[i-1] - ailen[i-1]; 697 rmax = PetscMax(rmax,ailen[i]); 698 if (fshift) { 699 ip = aj + ai[i] ; 700 ap = aa + ai[i] ; 701 N = ailen[i]; 702 for (j=0; j<N; j++) { 703 ip[j-fshift] = ip[j]; 704 ap[j-fshift] = ap[j]; 705 } 706 } 707 ai[i] = ai[i-1] + ailen[i-1]; 708 } 709 if (m) { 710 fshift += imax[m-1] - ailen[m-1]; 711 ai[m] = ai[m-1] + ailen[m-1]; 712 } 713 /* reset ilen and imax for each row */ 714 for (i=0; i<m; i++) { 715 ailen[i] = imax[i] = ai[i+1] - ai[i]; 716 } 717 a->nz = ai[m]; 718 if (fshift && a->nounused == -1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D, %D unneeded", m, A->cmap->n, fshift); 719 720 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 721 ierr = PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n,fshift,a->nz);CHKERRQ(ierr); 722 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 723 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 724 A->info.mallocs += a->reallocs; 725 a->reallocs = 0; 726 A->info.nz_unneeded = (double)fshift; 727 a->rmax = rmax; 728 729 /* check for zero rows. If found a large number of zero rows, use CompressedRow functions */ 730 ierr = Mat_CheckCompressedRow(A,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 731 A->same_nonzero = PETSC_TRUE; 732 733 ierr = MatAssemblyEnd_SeqAIJ_Inode(A,mode);CHKERRQ(ierr); 734 735 a->idiagvalid = PETSC_FALSE; 736 737 #if defined(PETSC_HAVE_CUDA) 738 a->GPUmatrix.resize(m,A->cmap->n,a->nz); 739 a->GPUmatrix.row_offsets.assign(a->i,a->i+m); 740 a->GPUmatrix.column_indices.assign(a->j,a->j+A->cmap->n); 741 a->GPUmatrix.values.assign(a->a,a->a+a->nz); 742 #endif 743 PetscFunctionReturn(0); 744 } 745 746 #undef __FUNCT__ 747 #define __FUNCT__ "MatRealPart_SeqAIJ" 748 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 749 { 750 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 751 PetscInt i,nz = a->nz; 752 MatScalar *aa = a->a; 753 754 PetscFunctionBegin; 755 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 756 PetscFunctionReturn(0); 757 } 758 759 #undef __FUNCT__ 760 #define __FUNCT__ "MatImaginaryPart_SeqAIJ" 761 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 762 { 763 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 764 PetscInt i,nz = a->nz; 765 MatScalar *aa = a->a; 766 767 PetscFunctionBegin; 768 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 769 PetscFunctionReturn(0); 770 } 771 772 #undef __FUNCT__ 773 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 774 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 775 { 776 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 777 PetscErrorCode ierr; 778 779 PetscFunctionBegin; 780 ierr = PetscMemzero(a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 781 PetscFunctionReturn(0); 782 } 783 784 #undef __FUNCT__ 785 #define __FUNCT__ "MatDestroy_SeqAIJ" 786 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 787 { 788 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 789 PetscErrorCode ierr; 790 791 PetscFunctionBegin; 792 #if defined(PETSC_USE_LOG) 793 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap->n,A->cmap->n,a->nz); 794 #endif 795 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 796 if (a->row) { 797 ierr = ISDestroy(a->row);CHKERRQ(ierr); 798 } 799 if (a->col) { 800 ierr = ISDestroy(a->col);CHKERRQ(ierr); 801 } 802 ierr = PetscFree(a->diag);CHKERRQ(ierr); 803 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 804 ierr = PetscFree3(a->idiag,a->mdiag,a->ssor_work);CHKERRQ(ierr); 805 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 806 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 807 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 808 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 809 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 810 if (a->XtoY) {ierr = MatDestroy(a->XtoY);CHKERRQ(ierr);} 811 if (a->compressedrow.checked && a->compressedrow.use){ierr = PetscFree2(a->compressedrow.i,a->compressedrow.rindex);} 812 813 ierr = MatDestroy_SeqAIJ_Inode(A);CHKERRQ(ierr); 814 815 ierr = PetscFree(a);CHKERRQ(ierr); 816 817 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 818 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 819 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 820 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 821 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqsbaij_C","",PETSC_NULL);CHKERRQ(ierr); 822 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 823 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqcsrperm_C","",PETSC_NULL);CHKERRQ(ierr); 824 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatIsTranspose_C","",PETSC_NULL);CHKERRQ(ierr); 825 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 826 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C","",PETSC_NULL);CHKERRQ(ierr); 827 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatReorderForNonzeroDiagonal_C","",PETSC_NULL);CHKERRQ(ierr); 828 PetscFunctionReturn(0); 829 } 830 831 #undef __FUNCT__ 832 #define __FUNCT__ "MatSetOption_SeqAIJ" 833 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscTruth flg) 834 { 835 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 836 PetscErrorCode ierr; 837 838 PetscFunctionBegin; 839 switch (op) { 840 case MAT_ROW_ORIENTED: 841 a->roworiented = flg; 842 break; 843 case MAT_KEEP_NONZERO_PATTERN: 844 a->keepnonzeropattern = flg; 845 break; 846 case MAT_NEW_NONZERO_LOCATIONS: 847 a->nonew = (flg ? 0 : 1); 848 break; 849 case MAT_NEW_NONZERO_LOCATION_ERR: 850 a->nonew = (flg ? -1 : 0); 851 break; 852 case MAT_NEW_NONZERO_ALLOCATION_ERR: 853 a->nonew = (flg ? -2 : 0); 854 break; 855 case MAT_UNUSED_NONZERO_LOCATION_ERR: 856 a->nounused = (flg ? -1 : 0); 857 break; 858 case MAT_IGNORE_ZERO_ENTRIES: 859 a->ignorezeroentries = flg; 860 break; 861 case MAT_USE_COMPRESSEDROW: 862 a->compressedrow.use = flg; 863 break; 864 case MAT_SYMMETRIC: 865 case MAT_STRUCTURALLY_SYMMETRIC: 866 case MAT_HERMITIAN: 867 case MAT_SYMMETRY_ETERNAL: 868 case MAT_NEW_DIAGONALS: 869 case MAT_IGNORE_OFF_PROC_ENTRIES: 870 case MAT_USE_HASH_TABLE: 871 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 872 break; 873 case MAT_USE_INODES: 874 /* Not an error because MatSetOption_SeqAIJ_Inode handles this one */ 875 break; 876 default: 877 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 878 } 879 ierr = MatSetOption_SeqAIJ_Inode(A,op,flg);CHKERRQ(ierr); 880 PetscFunctionReturn(0); 881 } 882 883 #undef __FUNCT__ 884 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 885 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 886 { 887 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 888 PetscErrorCode ierr; 889 PetscInt i,j,n,*ai=a->i,*aj=a->j,nz; 890 PetscScalar *aa=a->a,*x,zero=0.0; 891 892 PetscFunctionBegin; 893 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 894 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 895 896 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU){ 897 PetscInt *diag=a->diag; 898 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 899 for (i=0; i<n; i++) x[i] = aa[diag[i]]; 900 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 901 PetscFunctionReturn(0); 902 } 903 904 ierr = VecSet(v,zero);CHKERRQ(ierr); 905 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 906 for (i=0; i<n; i++) { 907 nz = ai[i+1] - ai[i]; 908 if (!nz) x[i] = 0.0; 909 for (j=ai[i]; j<ai[i+1]; j++){ 910 if (aj[j] == i) { 911 x[i] = aa[j]; 912 break; 913 } 914 } 915 } 916 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 917 PetscFunctionReturn(0); 918 } 919 920 #include "../src/mat/impls/aij/seq/ftn-kernels/fmult.h" 921 #undef __FUNCT__ 922 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 923 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 924 { 925 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 926 PetscScalar *x,*y; 927 PetscErrorCode ierr; 928 PetscInt m = A->rmap->n; 929 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 930 MatScalar *v; 931 PetscScalar alpha; 932 PetscInt n,i,j,*idx,*ii,*ridx=PETSC_NULL; 933 Mat_CompressedRow cprow = a->compressedrow; 934 PetscTruth usecprow = cprow.use; 935 #endif 936 937 PetscFunctionBegin; 938 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 939 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 940 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 941 942 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 943 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 944 #else 945 if (usecprow){ 946 m = cprow.nrows; 947 ii = cprow.i; 948 ridx = cprow.rindex; 949 } else { 950 ii = a->i; 951 } 952 for (i=0; i<m; i++) { 953 idx = a->j + ii[i] ; 954 v = a->a + ii[i] ; 955 n = ii[i+1] - ii[i]; 956 if (usecprow){ 957 alpha = x[ridx[i]]; 958 } else { 959 alpha = x[i]; 960 } 961 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 962 } 963 #endif 964 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 965 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 966 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 967 PetscFunctionReturn(0); 968 } 969 970 #undef __FUNCT__ 971 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 972 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 973 { 974 PetscErrorCode ierr; 975 976 PetscFunctionBegin; 977 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 978 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 979 PetscFunctionReturn(0); 980 } 981 982 #include "../src/mat/impls/aij/seq/ftn-kernels/fmult.h" 983 #undef __FUNCT__ 984 #define __FUNCT__ "MatMult_SeqAIJ" 985 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 986 { 987 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 988 PetscScalar *y; 989 const PetscScalar *x; 990 const MatScalar *aa; 991 PetscErrorCode ierr; 992 PetscInt m=A->rmap->n; 993 const PetscInt *aj,*ii,*ridx=PETSC_NULL; 994 PetscInt n,i,nonzerorow=0; 995 PetscScalar sum; 996 PetscTruth usecprow=a->compressedrow.use; 997 998 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 999 #pragma disjoint(*x,*y,*aa) 1000 #endif 1001 1002 PetscFunctionBegin; 1003 /*#if defined(PETSC_HAVE_CUDA) 1004 ierr = VecCUDACopyToGPU(xx);CHKERRQ(ierr); 1005 cusp::multiply(a->GPUmatrix,xx->GPUarray,yy->GPUarray); 1006 yy->valid_GPU_array = PETSC_CUDA_GPU; 1007 #else*/ 1008 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1009 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1010 aj = a->j; 1011 aa = a->a; 1012 ii = a->i; 1013 if (usecprow){ /* use compressed row format */ 1014 m = a->compressedrow.nrows; 1015 ii = a->compressedrow.i; 1016 ridx = a->compressedrow.rindex; 1017 for (i=0; i<m; i++){ 1018 n = ii[i+1] - ii[i]; 1019 aj = a->j + ii[i]; 1020 aa = a->a + ii[i]; 1021 sum = 0.0; 1022 nonzerorow += (n>0); 1023 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1024 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1025 y[*ridx++] = sum; 1026 } 1027 } else { /* do not use compressed row format */ 1028 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1029 fortranmultaij_(&m,x,ii,aj,aa,y); 1030 #else 1031 for (i=0; i<m; i++) { 1032 n = ii[i+1] - ii[i]; 1033 aj = a->j + ii[i]; 1034 aa = a->a + ii[i]; 1035 sum = 0.0; 1036 nonzerorow += (n>0); 1037 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1038 y[i] = sum; 1039 } 1040 #endif 1041 } 1042 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1043 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1044 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1045 /*#endif*/ 1046 PetscFunctionReturn(0); 1047 } 1048 1049 #include "../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h" 1050 #undef __FUNCT__ 1051 #define __FUNCT__ "MatMultAdd_SeqAIJ" 1052 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1053 { 1054 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1055 PetscScalar *x,*y,*z; 1056 const MatScalar *aa; 1057 PetscErrorCode ierr; 1058 PetscInt m = A->rmap->n,*aj,*ii; 1059 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1060 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 1061 PetscScalar sum; 1062 PetscTruth usecprow=a->compressedrow.use; 1063 #endif 1064 1065 PetscFunctionBegin; 1066 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1067 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1068 if (zz != yy) { 1069 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1070 } else { 1071 z = y; 1072 } 1073 1074 aj = a->j; 1075 aa = a->a; 1076 ii = a->i; 1077 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1078 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1079 #else 1080 if (usecprow){ /* use compressed row format */ 1081 if (zz != yy){ 1082 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1083 } 1084 m = a->compressedrow.nrows; 1085 ii = a->compressedrow.i; 1086 ridx = a->compressedrow.rindex; 1087 for (i=0; i<m; i++){ 1088 n = ii[i+1] - ii[i]; 1089 aj = a->j + ii[i]; 1090 aa = a->a + ii[i]; 1091 sum = y[*ridx]; 1092 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 1093 z[*ridx++] = sum; 1094 } 1095 } else { /* do not use compressed row format */ 1096 for (i=0; i<m; i++) { 1097 jrow = ii[i]; 1098 n = ii[i+1] - jrow; 1099 sum = y[i]; 1100 for (j=0; j<n; j++) { 1101 sum += aa[jrow]*x[aj[jrow]]; jrow++; 1102 } 1103 z[i] = sum; 1104 } 1105 } 1106 #endif 1107 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1108 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1109 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1110 if (zz != yy) { 1111 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1112 } 1113 PetscFunctionReturn(0); 1114 } 1115 1116 /* 1117 Adds diagonal pointers to sparse matrix structure. 1118 */ 1119 #undef __FUNCT__ 1120 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1121 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1122 { 1123 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1124 PetscErrorCode ierr; 1125 PetscInt i,j,m = A->rmap->n; 1126 1127 PetscFunctionBegin; 1128 if (!a->diag) { 1129 ierr = PetscMalloc(m*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 1130 ierr = PetscLogObjectMemory(A, m*sizeof(PetscInt));CHKERRQ(ierr); 1131 } 1132 for (i=0; i<A->rmap->n; i++) { 1133 a->diag[i] = a->i[i+1]; 1134 for (j=a->i[i]; j<a->i[i+1]; j++) { 1135 if (a->j[j] == i) { 1136 a->diag[i] = j; 1137 break; 1138 } 1139 } 1140 } 1141 PetscFunctionReturn(0); 1142 } 1143 1144 /* 1145 Checks for missing diagonals 1146 */ 1147 #undef __FUNCT__ 1148 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1149 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscTruth *missing,PetscInt *d) 1150 { 1151 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1152 PetscInt *diag,*jj = a->j,i; 1153 1154 PetscFunctionBegin; 1155 *missing = PETSC_FALSE; 1156 if (A->rmap->n > 0 && !jj) { 1157 *missing = PETSC_TRUE; 1158 if (d) *d = 0; 1159 PetscInfo(A,"Matrix has no entries therefor is missing diagonal"); 1160 } else { 1161 diag = a->diag; 1162 for (i=0; i<A->rmap->n; i++) { 1163 if (jj[diag[i]] != i) { 1164 *missing = PETSC_TRUE; 1165 if (d) *d = i; 1166 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1167 } 1168 } 1169 } 1170 PetscFunctionReturn(0); 1171 } 1172 1173 EXTERN_C_BEGIN 1174 #undef __FUNCT__ 1175 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1176 PetscErrorCode PETSCMAT_DLLEXPORT MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1177 { 1178 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1179 PetscErrorCode ierr; 1180 PetscInt i,*diag,m = A->rmap->n; 1181 MatScalar *v = a->a; 1182 PetscScalar *idiag,*mdiag; 1183 1184 PetscFunctionBegin; 1185 if (a->idiagvalid) PetscFunctionReturn(0); 1186 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1187 diag = a->diag; 1188 if (!a->idiag) { 1189 ierr = PetscMalloc3(m,PetscScalar,&a->idiag,m,PetscScalar,&a->mdiag,m,PetscScalar,&a->ssor_work);CHKERRQ(ierr); 1190 ierr = PetscLogObjectMemory(A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1191 v = a->a; 1192 } 1193 mdiag = a->mdiag; 1194 idiag = a->idiag; 1195 1196 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1197 for (i=0; i<m; i++) { 1198 mdiag[i] = v[diag[i]]; 1199 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1200 idiag[i] = 1.0/v[diag[i]]; 1201 } 1202 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1203 } else { 1204 for (i=0; i<m; i++) { 1205 mdiag[i] = v[diag[i]]; 1206 idiag[i] = omega/(fshift + v[diag[i]]); 1207 } 1208 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1209 } 1210 a->idiagvalid = PETSC_TRUE; 1211 PetscFunctionReturn(0); 1212 } 1213 EXTERN_C_END 1214 1215 #include "../src/mat/impls/aij/seq/ftn-kernels/frelax.h" 1216 #undef __FUNCT__ 1217 #define __FUNCT__ "MatSOR_SeqAIJ" 1218 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1219 { 1220 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1221 PetscScalar *x,d,sum,*t,scale; 1222 const MatScalar *v = a->a,*idiag=0,*mdiag; 1223 const PetscScalar *b, *bs,*xb, *ts; 1224 PetscErrorCode ierr; 1225 PetscInt n = A->cmap->n,m = A->rmap->n,i; 1226 const PetscInt *idx,*diag; 1227 1228 PetscFunctionBegin; 1229 its = its*lits; 1230 1231 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1232 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1233 a->fshift = fshift; 1234 a->omega = omega; 1235 1236 diag = a->diag; 1237 t = a->ssor_work; 1238 idiag = a->idiag; 1239 mdiag = a->mdiag; 1240 1241 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1242 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1243 CHKMEMQ; 1244 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1245 if (flag == SOR_APPLY_UPPER) { 1246 /* apply (U + D/omega) to the vector */ 1247 bs = b; 1248 for (i=0; i<m; i++) { 1249 d = fshift + mdiag[i]; 1250 n = a->i[i+1] - diag[i] - 1; 1251 idx = a->j + diag[i] + 1; 1252 v = a->a + diag[i] + 1; 1253 sum = b[i]*d/omega; 1254 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1255 x[i] = sum; 1256 } 1257 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1258 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1259 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1260 PetscFunctionReturn(0); 1261 } 1262 1263 if (flag == SOR_APPLY_LOWER) { 1264 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1265 } else if (flag & SOR_EISENSTAT) { 1266 /* Let A = L + U + D; where L is lower trianglar, 1267 U is upper triangular, E = D/omega; This routine applies 1268 1269 (L + E)^{-1} A (U + E)^{-1} 1270 1271 to a vector efficiently using Eisenstat's trick. 1272 */ 1273 scale = (2.0/omega) - 1.0; 1274 1275 /* x = (E + U)^{-1} b */ 1276 for (i=m-1; i>=0; i--) { 1277 n = a->i[i+1] - diag[i] - 1; 1278 idx = a->j + diag[i] + 1; 1279 v = a->a + diag[i] + 1; 1280 sum = b[i]; 1281 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1282 x[i] = sum*idiag[i]; 1283 } 1284 1285 /* t = b - (2*E - D)x */ 1286 v = a->a; 1287 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1288 1289 /* t = (E + L)^{-1}t */ 1290 ts = t; 1291 diag = a->diag; 1292 for (i=0; i<m; i++) { 1293 n = diag[i] - a->i[i]; 1294 idx = a->j + a->i[i]; 1295 v = a->a + a->i[i]; 1296 sum = t[i]; 1297 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1298 t[i] = sum*idiag[i]; 1299 /* x = x + t */ 1300 x[i] += t[i]; 1301 } 1302 1303 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1304 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1305 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1306 PetscFunctionReturn(0); 1307 } 1308 if (flag & SOR_ZERO_INITIAL_GUESS) { 1309 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1310 for (i=0; i<m; i++) { 1311 n = diag[i] - a->i[i]; 1312 idx = a->j + a->i[i]; 1313 v = a->a + a->i[i]; 1314 sum = b[i]; 1315 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1316 t[i] = sum; 1317 x[i] = sum*idiag[i]; 1318 } 1319 xb = t; 1320 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1321 } else xb = b; 1322 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1323 for (i=m-1; i>=0; i--) { 1324 n = a->i[i+1] - diag[i] - 1; 1325 idx = a->j + diag[i] + 1; 1326 v = a->a + diag[i] + 1; 1327 sum = xb[i]; 1328 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1329 if (xb == b) { 1330 x[i] = sum*idiag[i]; 1331 } else { 1332 x[i] = (1-omega)*x[i] + sum*idiag[i]; 1333 } 1334 } 1335 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1336 } 1337 its--; 1338 } 1339 while (its--) { 1340 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1341 for (i=0; i<m; i++) { 1342 n = a->i[i+1] - a->i[i]; 1343 idx = a->j + a->i[i]; 1344 v = a->a + a->i[i]; 1345 sum = b[i]; 1346 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1347 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1348 } 1349 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1350 } 1351 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1352 for (i=m-1; i>=0; i--) { 1353 n = a->i[i+1] - a->i[i]; 1354 idx = a->j + a->i[i]; 1355 v = a->a + a->i[i]; 1356 sum = b[i]; 1357 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1358 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1359 } 1360 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1361 } 1362 } 1363 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1364 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1365 CHKMEMQ; PetscFunctionReturn(0); 1366 } 1367 1368 1369 #undef __FUNCT__ 1370 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1371 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1372 { 1373 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1374 1375 PetscFunctionBegin; 1376 info->block_size = 1.0; 1377 info->nz_allocated = (double)a->maxnz; 1378 info->nz_used = (double)a->nz; 1379 info->nz_unneeded = (double)(a->maxnz - a->nz); 1380 info->assemblies = (double)A->num_ass; 1381 info->mallocs = (double)A->info.mallocs; 1382 info->memory = ((PetscObject)A)->mem; 1383 if (A->factortype) { 1384 info->fill_ratio_given = A->info.fill_ratio_given; 1385 info->fill_ratio_needed = A->info.fill_ratio_needed; 1386 info->factor_mallocs = A->info.factor_mallocs; 1387 } else { 1388 info->fill_ratio_given = 0; 1389 info->fill_ratio_needed = 0; 1390 info->factor_mallocs = 0; 1391 } 1392 PetscFunctionReturn(0); 1393 } 1394 1395 #undef __FUNCT__ 1396 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1397 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1398 { 1399 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1400 PetscInt i,m = A->rmap->n - 1,d = 0; 1401 PetscErrorCode ierr; 1402 PetscTruth missing; 1403 1404 PetscFunctionBegin; 1405 if (a->keepnonzeropattern) { 1406 for (i=0; i<N; i++) { 1407 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1408 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1409 } 1410 if (diag != 0.0) { 1411 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1412 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1413 for (i=0; i<N; i++) { 1414 a->a[a->diag[rows[i]]] = diag; 1415 } 1416 } 1417 A->same_nonzero = PETSC_TRUE; 1418 } else { 1419 if (diag != 0.0) { 1420 for (i=0; i<N; i++) { 1421 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1422 if (a->ilen[rows[i]] > 0) { 1423 a->ilen[rows[i]] = 1; 1424 a->a[a->i[rows[i]]] = diag; 1425 a->j[a->i[rows[i]]] = rows[i]; 1426 } else { /* in case row was completely empty */ 1427 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1428 } 1429 } 1430 } else { 1431 for (i=0; i<N; i++) { 1432 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1433 a->ilen[rows[i]] = 0; 1434 } 1435 } 1436 A->same_nonzero = PETSC_FALSE; 1437 } 1438 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1439 PetscFunctionReturn(0); 1440 } 1441 1442 #undef __FUNCT__ 1443 #define __FUNCT__ "MatGetRow_SeqAIJ" 1444 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1445 { 1446 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1447 PetscInt *itmp; 1448 1449 PetscFunctionBegin; 1450 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1451 1452 *nz = a->i[row+1] - a->i[row]; 1453 if (v) *v = a->a + a->i[row]; 1454 if (idx) { 1455 itmp = a->j + a->i[row]; 1456 if (*nz) { 1457 *idx = itmp; 1458 } 1459 else *idx = 0; 1460 } 1461 PetscFunctionReturn(0); 1462 } 1463 1464 /* remove this function? */ 1465 #undef __FUNCT__ 1466 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1467 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1468 { 1469 PetscFunctionBegin; 1470 PetscFunctionReturn(0); 1471 } 1472 1473 #undef __FUNCT__ 1474 #define __FUNCT__ "MatNorm_SeqAIJ" 1475 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1476 { 1477 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1478 MatScalar *v = a->a; 1479 PetscReal sum = 0.0; 1480 PetscErrorCode ierr; 1481 PetscInt i,j; 1482 1483 PetscFunctionBegin; 1484 if (type == NORM_FROBENIUS) { 1485 for (i=0; i<a->nz; i++) { 1486 #if defined(PETSC_USE_COMPLEX) 1487 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1488 #else 1489 sum += (*v)*(*v); v++; 1490 #endif 1491 } 1492 *nrm = sqrt(sum); 1493 } else if (type == NORM_1) { 1494 PetscReal *tmp; 1495 PetscInt *jj = a->j; 1496 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1497 ierr = PetscMemzero(tmp,A->cmap->n*sizeof(PetscReal));CHKERRQ(ierr); 1498 *nrm = 0.0; 1499 for (j=0; j<a->nz; j++) { 1500 tmp[*jj++] += PetscAbsScalar(*v); v++; 1501 } 1502 for (j=0; j<A->cmap->n; j++) { 1503 if (tmp[j] > *nrm) *nrm = tmp[j]; 1504 } 1505 ierr = PetscFree(tmp);CHKERRQ(ierr); 1506 } else if (type == NORM_INFINITY) { 1507 *nrm = 0.0; 1508 for (j=0; j<A->rmap->n; j++) { 1509 v = a->a + a->i[j]; 1510 sum = 0.0; 1511 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1512 sum += PetscAbsScalar(*v); v++; 1513 } 1514 if (sum > *nrm) *nrm = sum; 1515 } 1516 } else { 1517 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1518 } 1519 PetscFunctionReturn(0); 1520 } 1521 1522 #undef __FUNCT__ 1523 #define __FUNCT__ "MatTranspose_SeqAIJ" 1524 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 1525 { 1526 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1527 Mat C; 1528 PetscErrorCode ierr; 1529 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 1530 MatScalar *array = a->a; 1531 1532 PetscFunctionBegin; 1533 if (reuse == MAT_REUSE_MATRIX && A == *B && m != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1534 1535 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 1536 ierr = PetscMalloc((1+A->cmap->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1537 ierr = PetscMemzero(col,(1+A->cmap->n)*sizeof(PetscInt));CHKERRQ(ierr); 1538 1539 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1540 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1541 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 1542 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1543 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1544 ierr = PetscFree(col);CHKERRQ(ierr); 1545 } else { 1546 C = *B; 1547 } 1548 1549 for (i=0; i<m; i++) { 1550 len = ai[i+1]-ai[i]; 1551 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1552 array += len; 1553 aj += len; 1554 } 1555 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1556 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1557 1558 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 1559 *B = C; 1560 } else { 1561 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1562 } 1563 PetscFunctionReturn(0); 1564 } 1565 1566 EXTERN_C_BEGIN 1567 #undef __FUNCT__ 1568 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1569 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1570 { 1571 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1572 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1573 MatScalar *va,*vb; 1574 PetscErrorCode ierr; 1575 PetscInt ma,na,mb,nb, i; 1576 1577 PetscFunctionBegin; 1578 bij = (Mat_SeqAIJ *) B->data; 1579 1580 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1581 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1582 if (ma!=nb || na!=mb){ 1583 *f = PETSC_FALSE; 1584 PetscFunctionReturn(0); 1585 } 1586 aii = aij->i; bii = bij->i; 1587 adx = aij->j; bdx = bij->j; 1588 va = aij->a; vb = bij->a; 1589 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1590 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1591 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1592 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1593 1594 *f = PETSC_TRUE; 1595 for (i=0; i<ma; i++) { 1596 while (aptr[i]<aii[i+1]) { 1597 PetscInt idc,idr; 1598 PetscScalar vc,vr; 1599 /* column/row index/value */ 1600 idc = adx[aptr[i]]; 1601 idr = bdx[bptr[idc]]; 1602 vc = va[aptr[i]]; 1603 vr = vb[bptr[idc]]; 1604 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1605 *f = PETSC_FALSE; 1606 goto done; 1607 } else { 1608 aptr[i]++; 1609 if (B || i!=idc) bptr[idc]++; 1610 } 1611 } 1612 } 1613 done: 1614 ierr = PetscFree(aptr);CHKERRQ(ierr); 1615 if (B) { 1616 ierr = PetscFree(bptr);CHKERRQ(ierr); 1617 } 1618 PetscFunctionReturn(0); 1619 } 1620 EXTERN_C_END 1621 1622 EXTERN_C_BEGIN 1623 #undef __FUNCT__ 1624 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 1625 PetscErrorCode PETSCMAT_DLLEXPORT MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1626 { 1627 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1628 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1629 MatScalar *va,*vb; 1630 PetscErrorCode ierr; 1631 PetscInt ma,na,mb,nb, i; 1632 1633 PetscFunctionBegin; 1634 bij = (Mat_SeqAIJ *) B->data; 1635 1636 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1637 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1638 if (ma!=nb || na!=mb){ 1639 *f = PETSC_FALSE; 1640 PetscFunctionReturn(0); 1641 } 1642 aii = aij->i; bii = bij->i; 1643 adx = aij->j; bdx = bij->j; 1644 va = aij->a; vb = bij->a; 1645 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1646 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1647 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1648 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1649 1650 *f = PETSC_TRUE; 1651 for (i=0; i<ma; i++) { 1652 while (aptr[i]<aii[i+1]) { 1653 PetscInt idc,idr; 1654 PetscScalar vc,vr; 1655 /* column/row index/value */ 1656 idc = adx[aptr[i]]; 1657 idr = bdx[bptr[idc]]; 1658 vc = va[aptr[i]]; 1659 vr = vb[bptr[idc]]; 1660 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 1661 *f = PETSC_FALSE; 1662 goto done; 1663 } else { 1664 aptr[i]++; 1665 if (B || i!=idc) bptr[idc]++; 1666 } 1667 } 1668 } 1669 done: 1670 ierr = PetscFree(aptr);CHKERRQ(ierr); 1671 if (B) { 1672 ierr = PetscFree(bptr);CHKERRQ(ierr); 1673 } 1674 PetscFunctionReturn(0); 1675 } 1676 EXTERN_C_END 1677 1678 #undef __FUNCT__ 1679 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1680 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1681 { 1682 PetscErrorCode ierr; 1683 PetscFunctionBegin; 1684 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1685 PetscFunctionReturn(0); 1686 } 1687 1688 #undef __FUNCT__ 1689 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 1690 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1691 { 1692 PetscErrorCode ierr; 1693 PetscFunctionBegin; 1694 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1695 PetscFunctionReturn(0); 1696 } 1697 1698 #undef __FUNCT__ 1699 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1700 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1701 { 1702 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1703 PetscScalar *l,*r,x; 1704 MatScalar *v; 1705 PetscErrorCode ierr; 1706 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 1707 1708 PetscFunctionBegin; 1709 if (ll) { 1710 /* The local size is used so that VecMPI can be passed to this routine 1711 by MatDiagonalScale_MPIAIJ */ 1712 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1713 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1714 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1715 v = a->a; 1716 for (i=0; i<m; i++) { 1717 x = l[i]; 1718 M = a->i[i+1] - a->i[i]; 1719 for (j=0; j<M; j++) { (*v++) *= x;} 1720 } 1721 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1722 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1723 } 1724 if (rr) { 1725 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1726 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1727 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1728 v = a->a; jj = a->j; 1729 for (i=0; i<nz; i++) { 1730 (*v++) *= r[*jj++]; 1731 } 1732 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1733 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1734 } 1735 PetscFunctionReturn(0); 1736 } 1737 1738 #undef __FUNCT__ 1739 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1740 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1741 { 1742 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1743 PetscErrorCode ierr; 1744 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 1745 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1746 const PetscInt *irow,*icol; 1747 PetscInt nrows,ncols; 1748 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1749 MatScalar *a_new,*mat_a; 1750 Mat C; 1751 PetscTruth stride,sorted; 1752 1753 PetscFunctionBegin; 1754 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 1755 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1756 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 1757 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1758 1759 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1760 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1761 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1762 1763 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1764 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1765 if (stride && step == 1) { 1766 /* special case of contiguous rows */ 1767 ierr = PetscMalloc2(nrows,PetscInt,&lens,nrows,PetscInt,&starts);CHKERRQ(ierr); 1768 /* loop over new rows determining lens and starting points */ 1769 for (i=0; i<nrows; i++) { 1770 kstart = ai[irow[i]]; 1771 kend = kstart + ailen[irow[i]]; 1772 for (k=kstart; k<kend; k++) { 1773 if (aj[k] >= first) { 1774 starts[i] = k; 1775 break; 1776 } 1777 } 1778 sum = 0; 1779 while (k < kend) { 1780 if (aj[k++] >= first+ncols) break; 1781 sum++; 1782 } 1783 lens[i] = sum; 1784 } 1785 /* create submatrix */ 1786 if (scall == MAT_REUSE_MATRIX) { 1787 PetscInt n_cols,n_rows; 1788 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1789 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1790 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1791 C = *B; 1792 } else { 1793 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1794 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1795 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1796 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1797 } 1798 c = (Mat_SeqAIJ*)C->data; 1799 1800 /* loop over rows inserting into submatrix */ 1801 a_new = c->a; 1802 j_new = c->j; 1803 i_new = c->i; 1804 1805 for (i=0; i<nrows; i++) { 1806 ii = starts[i]; 1807 lensi = lens[i]; 1808 for (k=0; k<lensi; k++) { 1809 *j_new++ = aj[ii+k] - first; 1810 } 1811 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1812 a_new += lensi; 1813 i_new[i+1] = i_new[i] + lensi; 1814 c->ilen[i] = lensi; 1815 } 1816 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 1817 } else { 1818 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1819 ierr = PetscMalloc(oldcols*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1820 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1821 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1822 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1823 /* determine lens of each row */ 1824 for (i=0; i<nrows; i++) { 1825 kstart = ai[irow[i]]; 1826 kend = kstart + a->ilen[irow[i]]; 1827 lens[i] = 0; 1828 for (k=kstart; k<kend; k++) { 1829 if (smap[aj[k]]) { 1830 lens[i]++; 1831 } 1832 } 1833 } 1834 /* Create and fill new matrix */ 1835 if (scall == MAT_REUSE_MATRIX) { 1836 PetscTruth equal; 1837 1838 c = (Mat_SeqAIJ *)((*B)->data); 1839 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1840 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1841 if (!equal) { 1842 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1843 } 1844 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 1845 C = *B; 1846 } else { 1847 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1848 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1849 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1850 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1851 } 1852 c = (Mat_SeqAIJ *)(C->data); 1853 for (i=0; i<nrows; i++) { 1854 row = irow[i]; 1855 kstart = ai[row]; 1856 kend = kstart + a->ilen[row]; 1857 mat_i = c->i[i]; 1858 mat_j = c->j + mat_i; 1859 mat_a = c->a + mat_i; 1860 mat_ilen = c->ilen + i; 1861 for (k=kstart; k<kend; k++) { 1862 if ((tcol=smap[a->j[k]])) { 1863 *mat_j++ = tcol - 1; 1864 *mat_a++ = a->a[k]; 1865 (*mat_ilen)++; 1866 1867 } 1868 } 1869 } 1870 /* Free work space */ 1871 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1872 ierr = PetscFree(smap);CHKERRQ(ierr); 1873 ierr = PetscFree(lens);CHKERRQ(ierr); 1874 } 1875 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1876 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1877 1878 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1879 *B = C; 1880 PetscFunctionReturn(0); 1881 } 1882 1883 #undef __FUNCT__ 1884 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1885 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 1886 { 1887 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1888 PetscErrorCode ierr; 1889 Mat outA; 1890 PetscTruth row_identity,col_identity; 1891 1892 PetscFunctionBegin; 1893 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1894 1895 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1896 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1897 1898 outA = inA; 1899 outA->factortype = MAT_FACTOR_LU; 1900 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1901 if (a->row) { ierr = ISDestroy(a->row);CHKERRQ(ierr);} 1902 a->row = row; 1903 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1904 if (a->col) { ierr = ISDestroy(a->col);CHKERRQ(ierr);} 1905 a->col = col; 1906 1907 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1908 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1909 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1910 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1911 1912 if (!a->solve_work) { /* this matrix may have been factored before */ 1913 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1914 ierr = PetscLogObjectMemory(inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 1915 } 1916 1917 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1918 if (row_identity && col_identity) { 1919 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 1920 } else { 1921 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 1922 } 1923 PetscFunctionReturn(0); 1924 } 1925 1926 #undef __FUNCT__ 1927 #define __FUNCT__ "MatScale_SeqAIJ" 1928 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1929 { 1930 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1931 PetscScalar oalpha = alpha; 1932 PetscErrorCode ierr; 1933 PetscBLASInt one = 1,bnz = PetscBLASIntCast(a->nz); 1934 1935 PetscFunctionBegin; 1936 BLASscal_(&bnz,&oalpha,a->a,&one); 1937 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1938 PetscFunctionReturn(0); 1939 } 1940 1941 #undef __FUNCT__ 1942 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1943 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1944 { 1945 PetscErrorCode ierr; 1946 PetscInt i; 1947 1948 PetscFunctionBegin; 1949 if (scall == MAT_INITIAL_MATRIX) { 1950 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1951 } 1952 1953 for (i=0; i<n; i++) { 1954 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1955 } 1956 PetscFunctionReturn(0); 1957 } 1958 1959 #undef __FUNCT__ 1960 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1961 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1962 { 1963 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1964 PetscErrorCode ierr; 1965 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 1966 const PetscInt *idx; 1967 PetscInt start,end,*ai,*aj; 1968 PetscBT table; 1969 1970 PetscFunctionBegin; 1971 m = A->rmap->n; 1972 ai = a->i; 1973 aj = a->j; 1974 1975 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1976 1977 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1978 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1979 1980 for (i=0; i<is_max; i++) { 1981 /* Initialize the two local arrays */ 1982 isz = 0; 1983 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1984 1985 /* Extract the indices, assume there can be duplicate entries */ 1986 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1987 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1988 1989 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1990 for (j=0; j<n ; ++j){ 1991 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1992 } 1993 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1994 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1995 1996 k = 0; 1997 for (j=0; j<ov; j++){ /* for each overlap */ 1998 n = isz; 1999 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 2000 row = nidx[k]; 2001 start = ai[row]; 2002 end = ai[row+1]; 2003 for (l = start; l<end ; l++){ 2004 val = aj[l] ; 2005 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 2006 } 2007 } 2008 } 2009 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 2010 } 2011 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 2012 ierr = PetscFree(nidx);CHKERRQ(ierr); 2013 PetscFunctionReturn(0); 2014 } 2015 2016 /* -------------------------------------------------------------- */ 2017 #undef __FUNCT__ 2018 #define __FUNCT__ "MatPermute_SeqAIJ" 2019 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2020 { 2021 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2022 PetscErrorCode ierr; 2023 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2024 const PetscInt *row,*col; 2025 PetscInt *cnew,j,*lens; 2026 IS icolp,irowp; 2027 PetscInt *cwork = PETSC_NULL; 2028 PetscScalar *vwork = PETSC_NULL; 2029 2030 PetscFunctionBegin; 2031 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2032 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2033 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2034 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2035 2036 /* determine lengths of permuted rows */ 2037 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2038 for (i=0; i<m; i++) { 2039 lens[row[i]] = a->i[i+1] - a->i[i]; 2040 } 2041 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 2042 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2043 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2044 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2045 ierr = PetscFree(lens);CHKERRQ(ierr); 2046 2047 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2048 for (i=0; i<m; i++) { 2049 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2050 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 2051 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2052 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2053 } 2054 ierr = PetscFree(cnew);CHKERRQ(ierr); 2055 (*B)->assembled = PETSC_FALSE; 2056 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2057 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2058 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2059 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2060 ierr = ISDestroy(irowp);CHKERRQ(ierr); 2061 ierr = ISDestroy(icolp);CHKERRQ(ierr); 2062 PetscFunctionReturn(0); 2063 } 2064 2065 #undef __FUNCT__ 2066 #define __FUNCT__ "MatCopy_SeqAIJ" 2067 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2068 { 2069 PetscErrorCode ierr; 2070 2071 PetscFunctionBegin; 2072 /* If the two matrices have the same copy implementation, use fast copy. */ 2073 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2074 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2075 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2076 2077 if (a->i[A->rmap->n] != b->i[B->rmap->n]) { 2078 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 2079 } 2080 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2081 } else { 2082 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2083 } 2084 PetscFunctionReturn(0); 2085 } 2086 2087 #undef __FUNCT__ 2088 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 2089 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 2090 { 2091 PetscErrorCode ierr; 2092 2093 PetscFunctionBegin; 2094 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2095 PetscFunctionReturn(0); 2096 } 2097 2098 #undef __FUNCT__ 2099 #define __FUNCT__ "MatGetArray_SeqAIJ" 2100 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2101 { 2102 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2103 PetscFunctionBegin; 2104 *array = a->a; 2105 PetscFunctionReturn(0); 2106 } 2107 2108 #undef __FUNCT__ 2109 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2110 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2111 { 2112 PetscFunctionBegin; 2113 PetscFunctionReturn(0); 2114 } 2115 2116 #undef __FUNCT__ 2117 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2118 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2119 { 2120 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2121 PetscErrorCode ierr; 2122 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2123 PetscScalar dx,*y,*xx,*w3_array; 2124 PetscScalar *vscale_array; 2125 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2126 Vec w1,w2,w3; 2127 void *fctx = coloring->fctx; 2128 PetscTruth flg = PETSC_FALSE; 2129 2130 PetscFunctionBegin; 2131 if (!coloring->w1) { 2132 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2133 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2134 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2135 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2136 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2137 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2138 } 2139 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2140 2141 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2142 ierr = PetscOptionsGetTruth(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,PETSC_NULL);CHKERRQ(ierr); 2143 if (flg) { 2144 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2145 } else { 2146 PetscTruth assembled; 2147 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2148 if (assembled) { 2149 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2150 } 2151 } 2152 2153 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2154 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2155 2156 /* 2157 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2158 coloring->F for the coarser grids from the finest 2159 */ 2160 if (coloring->F) { 2161 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2162 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2163 if (m1 != m2) { 2164 coloring->F = 0; 2165 } 2166 } 2167 2168 if (coloring->F) { 2169 w1 = coloring->F; 2170 coloring->F = 0; 2171 } else { 2172 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2173 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2174 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2175 } 2176 2177 /* 2178 Compute all the scale factors and share with other processors 2179 */ 2180 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2181 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2182 for (k=0; k<coloring->ncolors; k++) { 2183 /* 2184 Loop over each column associated with color adding the 2185 perturbation to the vector w3. 2186 */ 2187 for (l=0; l<coloring->ncolumns[k]; l++) { 2188 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2189 dx = xx[col]; 2190 if (dx == 0.0) dx = 1.0; 2191 #if !defined(PETSC_USE_COMPLEX) 2192 if (dx < umin && dx >= 0.0) dx = umin; 2193 else if (dx < 0.0 && dx > -umin) dx = -umin; 2194 #else 2195 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2196 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2197 #endif 2198 dx *= epsilon; 2199 vscale_array[col] = 1.0/dx; 2200 } 2201 } 2202 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2203 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2204 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2205 2206 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2207 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2208 2209 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2210 else vscaleforrow = coloring->columnsforrow; 2211 2212 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2213 /* 2214 Loop over each color 2215 */ 2216 for (k=0; k<coloring->ncolors; k++) { 2217 coloring->currentcolor = k; 2218 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2219 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2220 /* 2221 Loop over each column associated with color adding the 2222 perturbation to the vector w3. 2223 */ 2224 for (l=0; l<coloring->ncolumns[k]; l++) { 2225 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2226 dx = xx[col]; 2227 if (dx == 0.0) dx = 1.0; 2228 #if !defined(PETSC_USE_COMPLEX) 2229 if (dx < umin && dx >= 0.0) dx = umin; 2230 else if (dx < 0.0 && dx > -umin) dx = -umin; 2231 #else 2232 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2233 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2234 #endif 2235 dx *= epsilon; 2236 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2237 w3_array[col] += dx; 2238 } 2239 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2240 2241 /* 2242 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2243 */ 2244 2245 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2246 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2247 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2248 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2249 2250 /* 2251 Loop over rows of vector, putting results into Jacobian matrix 2252 */ 2253 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2254 for (l=0; l<coloring->nrows[k]; l++) { 2255 row = coloring->rows[k][l]; 2256 col = coloring->columnsforrow[k][l]; 2257 y[row] *= vscale_array[vscaleforrow[k][l]]; 2258 srow = row + start; 2259 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2260 } 2261 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2262 } 2263 coloring->currentcolor = k; 2264 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2265 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2266 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2267 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2268 PetscFunctionReturn(0); 2269 } 2270 2271 #undef __FUNCT__ 2272 #define __FUNCT__ "MatAXPY_SeqAIJ" 2273 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2274 { 2275 PetscErrorCode ierr; 2276 PetscInt i; 2277 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2278 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2279 2280 PetscFunctionBegin; 2281 if (str == SAME_NONZERO_PATTERN) { 2282 PetscScalar alpha = a; 2283 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2284 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2285 if (y->xtoy && y->XtoY != X) { 2286 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2287 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2288 } 2289 if (!y->xtoy) { /* get xtoy */ 2290 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2291 y->XtoY = X; 2292 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2293 } 2294 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2295 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2296 } else { 2297 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2298 } 2299 PetscFunctionReturn(0); 2300 } 2301 2302 #undef __FUNCT__ 2303 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2304 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2305 { 2306 PetscErrorCode ierr; 2307 2308 PetscFunctionBegin; 2309 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 2310 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 2311 PetscFunctionReturn(0); 2312 } 2313 2314 #undef __FUNCT__ 2315 #define __FUNCT__ "MatConjugate_SeqAIJ" 2316 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2317 { 2318 #if defined(PETSC_USE_COMPLEX) 2319 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2320 PetscInt i,nz; 2321 PetscScalar *a; 2322 2323 PetscFunctionBegin; 2324 nz = aij->nz; 2325 a = aij->a; 2326 for (i=0; i<nz; i++) { 2327 a[i] = PetscConj(a[i]); 2328 } 2329 #else 2330 PetscFunctionBegin; 2331 #endif 2332 PetscFunctionReturn(0); 2333 } 2334 2335 #undef __FUNCT__ 2336 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2337 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2338 { 2339 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2340 PetscErrorCode ierr; 2341 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2342 PetscReal atmp; 2343 PetscScalar *x; 2344 MatScalar *aa; 2345 2346 PetscFunctionBegin; 2347 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2348 aa = a->a; 2349 ai = a->i; 2350 aj = a->j; 2351 2352 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2353 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2354 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2355 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2356 for (i=0; i<m; i++) { 2357 ncols = ai[1] - ai[0]; ai++; 2358 x[i] = 0.0; 2359 for (j=0; j<ncols; j++){ 2360 atmp = PetscAbsScalar(*aa); 2361 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2362 aa++; aj++; 2363 } 2364 } 2365 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2366 PetscFunctionReturn(0); 2367 } 2368 2369 #undef __FUNCT__ 2370 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2371 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2372 { 2373 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2374 PetscErrorCode ierr; 2375 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2376 PetscScalar *x; 2377 MatScalar *aa; 2378 2379 PetscFunctionBegin; 2380 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2381 aa = a->a; 2382 ai = a->i; 2383 aj = a->j; 2384 2385 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2386 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2387 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2388 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2389 for (i=0; i<m; i++) { 2390 ncols = ai[1] - ai[0]; ai++; 2391 if (ncols == A->cmap->n) { /* row is dense */ 2392 x[i] = *aa; if (idx) idx[i] = 0; 2393 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2394 x[i] = 0.0; 2395 if (idx) { 2396 idx[i] = 0; /* in case ncols is zero */ 2397 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2398 if (aj[j] > j) { 2399 idx[i] = j; 2400 break; 2401 } 2402 } 2403 } 2404 } 2405 for (j=0; j<ncols; j++){ 2406 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2407 aa++; aj++; 2408 } 2409 } 2410 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2411 PetscFunctionReturn(0); 2412 } 2413 2414 #undef __FUNCT__ 2415 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2416 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2417 { 2418 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2419 PetscErrorCode ierr; 2420 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2421 PetscReal atmp; 2422 PetscScalar *x; 2423 MatScalar *aa; 2424 2425 PetscFunctionBegin; 2426 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2427 aa = a->a; 2428 ai = a->i; 2429 aj = a->j; 2430 2431 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2432 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2433 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2434 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2435 for (i=0; i<m; i++) { 2436 ncols = ai[1] - ai[0]; ai++; 2437 if (ncols) { 2438 /* Get first nonzero */ 2439 for(j = 0; j < ncols; j++) { 2440 atmp = PetscAbsScalar(aa[j]); 2441 if (atmp > 1.0e-12) {x[i] = atmp; if (idx) idx[i] = aj[j]; break;} 2442 } 2443 if (j == ncols) {x[i] = *aa; if (idx) idx[i] = *aj;} 2444 } else { 2445 x[i] = 0.0; if (idx) idx[i] = 0; 2446 } 2447 for(j = 0; j < ncols; j++) { 2448 atmp = PetscAbsScalar(*aa); 2449 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2450 aa++; aj++; 2451 } 2452 } 2453 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2454 PetscFunctionReturn(0); 2455 } 2456 2457 #undef __FUNCT__ 2458 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2459 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2460 { 2461 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2462 PetscErrorCode ierr; 2463 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2464 PetscScalar *x; 2465 MatScalar *aa; 2466 2467 PetscFunctionBegin; 2468 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2469 aa = a->a; 2470 ai = a->i; 2471 aj = a->j; 2472 2473 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2474 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2475 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2476 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2477 for (i=0; i<m; i++) { 2478 ncols = ai[1] - ai[0]; ai++; 2479 if (ncols == A->cmap->n) { /* row is dense */ 2480 x[i] = *aa; if (idx) idx[i] = 0; 2481 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2482 x[i] = 0.0; 2483 if (idx) { /* find first implicit 0.0 in the row */ 2484 idx[i] = 0; /* in case ncols is zero */ 2485 for (j=0;j<ncols;j++) { 2486 if (aj[j] > j) { 2487 idx[i] = j; 2488 break; 2489 } 2490 } 2491 } 2492 } 2493 for (j=0; j<ncols; j++){ 2494 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2495 aa++; aj++; 2496 } 2497 } 2498 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2499 PetscFunctionReturn(0); 2500 } 2501 extern PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 2502 /* -------------------------------------------------------------------*/ 2503 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2504 MatGetRow_SeqAIJ, 2505 MatRestoreRow_SeqAIJ, 2506 MatMult_SeqAIJ, 2507 /* 4*/ MatMultAdd_SeqAIJ, 2508 MatMultTranspose_SeqAIJ, 2509 MatMultTransposeAdd_SeqAIJ, 2510 0, 2511 0, 2512 0, 2513 /*10*/ 0, 2514 MatLUFactor_SeqAIJ, 2515 0, 2516 MatSOR_SeqAIJ, 2517 MatTranspose_SeqAIJ, 2518 /*15*/ MatGetInfo_SeqAIJ, 2519 MatEqual_SeqAIJ, 2520 MatGetDiagonal_SeqAIJ, 2521 MatDiagonalScale_SeqAIJ, 2522 MatNorm_SeqAIJ, 2523 /*20*/ 0, 2524 MatAssemblyEnd_SeqAIJ, 2525 MatSetOption_SeqAIJ, 2526 MatZeroEntries_SeqAIJ, 2527 /*24*/ MatZeroRows_SeqAIJ, 2528 0, 2529 0, 2530 0, 2531 0, 2532 /*29*/ MatSetUpPreallocation_SeqAIJ, 2533 0, 2534 0, 2535 MatGetArray_SeqAIJ, 2536 MatRestoreArray_SeqAIJ, 2537 /*34*/ MatDuplicate_SeqAIJ, 2538 0, 2539 0, 2540 MatILUFactor_SeqAIJ, 2541 0, 2542 /*39*/ MatAXPY_SeqAIJ, 2543 MatGetSubMatrices_SeqAIJ, 2544 MatIncreaseOverlap_SeqAIJ, 2545 MatGetValues_SeqAIJ, 2546 MatCopy_SeqAIJ, 2547 /*44*/ MatGetRowMax_SeqAIJ, 2548 MatScale_SeqAIJ, 2549 0, 2550 MatDiagonalSet_SeqAIJ, 2551 0, 2552 /*49*/ MatSetBlockSize_SeqAIJ, 2553 MatGetRowIJ_SeqAIJ, 2554 MatRestoreRowIJ_SeqAIJ, 2555 MatGetColumnIJ_SeqAIJ, 2556 MatRestoreColumnIJ_SeqAIJ, 2557 /*54*/ MatFDColoringCreate_SeqAIJ, 2558 0, 2559 0, 2560 MatPermute_SeqAIJ, 2561 0, 2562 /*59*/ 0, 2563 MatDestroy_SeqAIJ, 2564 MatView_SeqAIJ, 2565 0, 2566 0, 2567 /*64*/ 0, 2568 0, 2569 0, 2570 0, 2571 0, 2572 /*69*/ MatGetRowMaxAbs_SeqAIJ, 2573 MatGetRowMinAbs_SeqAIJ, 2574 0, 2575 MatSetColoring_SeqAIJ, 2576 #if defined(PETSC_HAVE_ADIC) 2577 MatSetValuesAdic_SeqAIJ, 2578 #else 2579 0, 2580 #endif 2581 /*74*/ MatSetValuesAdifor_SeqAIJ, 2582 MatFDColoringApply_AIJ, 2583 0, 2584 0, 2585 0, 2586 /*79*/ 0, 2587 0, 2588 0, 2589 0, 2590 MatLoad_SeqAIJ, 2591 /*84*/ MatIsSymmetric_SeqAIJ, 2592 MatIsHermitian_SeqAIJ, 2593 0, 2594 0, 2595 0, 2596 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 2597 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2598 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2599 MatPtAP_Basic, 2600 MatPtAPSymbolic_SeqAIJ, 2601 /*94*/ MatPtAPNumeric_SeqAIJ, 2602 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2603 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2604 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2605 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2606 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 2607 0, 2608 0, 2609 MatConjugate_SeqAIJ, 2610 0, 2611 /*104*/MatSetValuesRow_SeqAIJ, 2612 MatRealPart_SeqAIJ, 2613 MatImaginaryPart_SeqAIJ, 2614 0, 2615 0, 2616 /*109*/0, 2617 0, 2618 MatGetRowMin_SeqAIJ, 2619 0, 2620 MatMissingDiagonal_SeqAIJ, 2621 /*114*/0, 2622 0, 2623 0, 2624 0, 2625 0, 2626 /*119*/0, 2627 0, 2628 0, 2629 0, 2630 MatLoadnew_SeqAIJ 2631 }; 2632 2633 EXTERN_C_BEGIN 2634 #undef __FUNCT__ 2635 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2636 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2637 { 2638 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2639 PetscInt i,nz,n; 2640 2641 PetscFunctionBegin; 2642 2643 nz = aij->maxnz; 2644 n = mat->rmap->n; 2645 for (i=0; i<nz; i++) { 2646 aij->j[i] = indices[i]; 2647 } 2648 aij->nz = nz; 2649 for (i=0; i<n; i++) { 2650 aij->ilen[i] = aij->imax[i]; 2651 } 2652 2653 PetscFunctionReturn(0); 2654 } 2655 EXTERN_C_END 2656 2657 #undef __FUNCT__ 2658 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2659 /*@ 2660 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2661 in the matrix. 2662 2663 Input Parameters: 2664 + mat - the SeqAIJ matrix 2665 - indices - the column indices 2666 2667 Level: advanced 2668 2669 Notes: 2670 This can be called if you have precomputed the nonzero structure of the 2671 matrix and want to provide it to the matrix object to improve the performance 2672 of the MatSetValues() operation. 2673 2674 You MUST have set the correct numbers of nonzeros per row in the call to 2675 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2676 2677 MUST be called before any calls to MatSetValues(); 2678 2679 The indices should start with zero, not one. 2680 2681 @*/ 2682 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2683 { 2684 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2685 2686 PetscFunctionBegin; 2687 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2688 PetscValidPointer(indices,2); 2689 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2690 if (f) { 2691 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2692 } else { 2693 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2694 } 2695 PetscFunctionReturn(0); 2696 } 2697 2698 /* ----------------------------------------------------------------------------------------*/ 2699 2700 EXTERN_C_BEGIN 2701 #undef __FUNCT__ 2702 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2703 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2704 { 2705 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2706 PetscErrorCode ierr; 2707 size_t nz = aij->i[mat->rmap->n]; 2708 2709 PetscFunctionBegin; 2710 if (aij->nonew != 1) { 2711 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2712 } 2713 2714 /* allocate space for values if not already there */ 2715 if (!aij->saved_values) { 2716 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2717 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2718 } 2719 2720 /* copy values over */ 2721 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2722 PetscFunctionReturn(0); 2723 } 2724 EXTERN_C_END 2725 2726 #undef __FUNCT__ 2727 #define __FUNCT__ "MatStoreValues" 2728 /*@ 2729 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2730 example, reuse of the linear part of a Jacobian, while recomputing the 2731 nonlinear portion. 2732 2733 Collect on Mat 2734 2735 Input Parameters: 2736 . mat - the matrix (currently only AIJ matrices support this option) 2737 2738 Level: advanced 2739 2740 Common Usage, with SNESSolve(): 2741 $ Create Jacobian matrix 2742 $ Set linear terms into matrix 2743 $ Apply boundary conditions to matrix, at this time matrix must have 2744 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2745 $ boundary conditions again will not change the nonzero structure 2746 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2747 $ ierr = MatStoreValues(mat); 2748 $ Call SNESSetJacobian() with matrix 2749 $ In your Jacobian routine 2750 $ ierr = MatRetrieveValues(mat); 2751 $ Set nonlinear terms in matrix 2752 2753 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2754 $ // build linear portion of Jacobian 2755 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2756 $ ierr = MatStoreValues(mat); 2757 $ loop over nonlinear iterations 2758 $ ierr = MatRetrieveValues(mat); 2759 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2760 $ // call MatAssemblyBegin/End() on matrix 2761 $ Solve linear system with Jacobian 2762 $ endloop 2763 2764 Notes: 2765 Matrix must already be assemblied before calling this routine 2766 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2767 calling this routine. 2768 2769 When this is called multiple times it overwrites the previous set of stored values 2770 and does not allocated additional space. 2771 2772 .seealso: MatRetrieveValues() 2773 2774 @*/ 2775 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2776 { 2777 PetscErrorCode ierr,(*f)(Mat); 2778 2779 PetscFunctionBegin; 2780 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2781 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2782 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2783 2784 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2785 if (f) { 2786 ierr = (*f)(mat);CHKERRQ(ierr); 2787 } else { 2788 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2789 } 2790 PetscFunctionReturn(0); 2791 } 2792 2793 EXTERN_C_BEGIN 2794 #undef __FUNCT__ 2795 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2796 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2797 { 2798 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2799 PetscErrorCode ierr; 2800 PetscInt nz = aij->i[mat->rmap->n]; 2801 2802 PetscFunctionBegin; 2803 if (aij->nonew != 1) { 2804 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2805 } 2806 if (!aij->saved_values) { 2807 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2808 } 2809 /* copy values over */ 2810 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2811 PetscFunctionReturn(0); 2812 } 2813 EXTERN_C_END 2814 2815 #undef __FUNCT__ 2816 #define __FUNCT__ "MatRetrieveValues" 2817 /*@ 2818 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2819 example, reuse of the linear part of a Jacobian, while recomputing the 2820 nonlinear portion. 2821 2822 Collect on Mat 2823 2824 Input Parameters: 2825 . mat - the matrix (currently on AIJ matrices support this option) 2826 2827 Level: advanced 2828 2829 .seealso: MatStoreValues() 2830 2831 @*/ 2832 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2833 { 2834 PetscErrorCode ierr,(*f)(Mat); 2835 2836 PetscFunctionBegin; 2837 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2838 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2839 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2840 2841 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2842 if (f) { 2843 ierr = (*f)(mat);CHKERRQ(ierr); 2844 } else { 2845 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2846 } 2847 PetscFunctionReturn(0); 2848 } 2849 2850 2851 /* --------------------------------------------------------------------------------*/ 2852 #undef __FUNCT__ 2853 #define __FUNCT__ "MatCreateSeqAIJ" 2854 /*@C 2855 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2856 (the default parallel PETSc format). For good matrix assembly performance 2857 the user should preallocate the matrix storage by setting the parameter nz 2858 (or the array nnz). By setting these parameters accurately, performance 2859 during matrix assembly can be increased by more than a factor of 50. 2860 2861 Collective on MPI_Comm 2862 2863 Input Parameters: 2864 + comm - MPI communicator, set to PETSC_COMM_SELF 2865 . m - number of rows 2866 . n - number of columns 2867 . nz - number of nonzeros per row (same for all rows) 2868 - nnz - array containing the number of nonzeros in the various rows 2869 (possibly different for each row) or PETSC_NULL 2870 2871 Output Parameter: 2872 . A - the matrix 2873 2874 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2875 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2876 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2877 2878 Notes: 2879 If nnz is given then nz is ignored 2880 2881 The AIJ format (also called the Yale sparse matrix format or 2882 compressed row storage), is fully compatible with standard Fortran 77 2883 storage. That is, the stored row and column indices can begin at 2884 either one (as in Fortran) or zero. See the users' manual for details. 2885 2886 Specify the preallocated storage with either nz or nnz (not both). 2887 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2888 allocation. For large problems you MUST preallocate memory or you 2889 will get TERRIBLE performance, see the users' manual chapter on matrices. 2890 2891 By default, this format uses inodes (identical nodes) when possible, to 2892 improve numerical efficiency of matrix-vector products and solves. We 2893 search for consecutive rows with the same nonzero structure, thereby 2894 reusing matrix information to achieve increased efficiency. 2895 2896 Options Database Keys: 2897 + -mat_no_inode - Do not use inodes 2898 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2899 - -mat_aij_oneindex - Internally use indexing starting at 1 2900 rather than 0. Note that when calling MatSetValues(), 2901 the user still MUST index entries starting at 0! 2902 2903 Level: intermediate 2904 2905 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2906 2907 @*/ 2908 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2909 { 2910 PetscErrorCode ierr; 2911 2912 PetscFunctionBegin; 2913 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2914 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2915 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2916 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2917 PetscFunctionReturn(0); 2918 } 2919 2920 #undef __FUNCT__ 2921 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2922 /*@C 2923 MatSeqAIJSetPreallocation - For good matrix assembly performance 2924 the user should preallocate the matrix storage by setting the parameter nz 2925 (or the array nnz). By setting these parameters accurately, performance 2926 during matrix assembly can be increased by more than a factor of 50. 2927 2928 Collective on MPI_Comm 2929 2930 Input Parameters: 2931 + B - The matrix-free 2932 . nz - number of nonzeros per row (same for all rows) 2933 - nnz - array containing the number of nonzeros in the various rows 2934 (possibly different for each row) or PETSC_NULL 2935 2936 Notes: 2937 If nnz is given then nz is ignored 2938 2939 The AIJ format (also called the Yale sparse matrix format or 2940 compressed row storage), is fully compatible with standard Fortran 77 2941 storage. That is, the stored row and column indices can begin at 2942 either one (as in Fortran) or zero. See the users' manual for details. 2943 2944 Specify the preallocated storage with either nz or nnz (not both). 2945 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2946 allocation. For large problems you MUST preallocate memory or you 2947 will get TERRIBLE performance, see the users' manual chapter on matrices. 2948 2949 You can call MatGetInfo() to get information on how effective the preallocation was; 2950 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2951 You can also run with the option -info and look for messages with the string 2952 malloc in them to see if additional memory allocation was needed. 2953 2954 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2955 entries or columns indices 2956 2957 By default, this format uses inodes (identical nodes) when possible, to 2958 improve numerical efficiency of matrix-vector products and solves. We 2959 search for consecutive rows with the same nonzero structure, thereby 2960 reusing matrix information to achieve increased efficiency. 2961 2962 Options Database Keys: 2963 + -mat_no_inode - Do not use inodes 2964 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2965 - -mat_aij_oneindex - Internally use indexing starting at 1 2966 rather than 0. Note that when calling MatSetValues(), 2967 the user still MUST index entries starting at 0! 2968 2969 Level: intermediate 2970 2971 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 2972 2973 @*/ 2974 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2975 { 2976 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2977 2978 PetscFunctionBegin; 2979 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2980 if (f) { 2981 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2982 } 2983 PetscFunctionReturn(0); 2984 } 2985 2986 EXTERN_C_BEGIN 2987 #undef __FUNCT__ 2988 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2989 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2990 { 2991 Mat_SeqAIJ *b; 2992 PetscTruth skipallocation = PETSC_FALSE; 2993 PetscErrorCode ierr; 2994 PetscInt i; 2995 2996 PetscFunctionBegin; 2997 2998 if (nz == MAT_SKIP_ALLOCATION) { 2999 skipallocation = PETSC_TRUE; 3000 nz = 0; 3001 } 3002 3003 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 3004 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 3005 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3006 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3007 3008 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3009 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3010 if (nnz) { 3011 for (i=0; i<B->rmap->n; i++) { 3012 if (nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 3013 if (nnz[i] > B->cmap->n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->cmap->n); 3014 } 3015 } 3016 3017 B->preallocated = PETSC_TRUE; 3018 b = (Mat_SeqAIJ*)B->data; 3019 3020 if (!skipallocation) { 3021 if (!b->imax) { 3022 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3023 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3024 } 3025 if (!nnz) { 3026 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3027 else if (nz <= 0) nz = 1; 3028 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3029 nz = nz*B->rmap->n; 3030 } else { 3031 nz = 0; 3032 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3033 } 3034 /* b->ilen will count nonzeros in each row so far. */ 3035 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3036 3037 /* allocate the matrix space */ 3038 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3039 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3040 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3041 b->i[0] = 0; 3042 for (i=1; i<B->rmap->n+1; i++) { 3043 b->i[i] = b->i[i-1] + b->imax[i-1]; 3044 } 3045 b->singlemalloc = PETSC_TRUE; 3046 b->free_a = PETSC_TRUE; 3047 b->free_ij = PETSC_TRUE; 3048 } else { 3049 b->free_a = PETSC_FALSE; 3050 b->free_ij = PETSC_FALSE; 3051 } 3052 3053 b->nz = 0; 3054 b->maxnz = nz; 3055 B->info.nz_unneeded = (double)b->maxnz; 3056 PetscFunctionReturn(0); 3057 } 3058 EXTERN_C_END 3059 3060 #undef __FUNCT__ 3061 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3062 /*@ 3063 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3064 3065 Input Parameters: 3066 + B - the matrix 3067 . i - the indices into j for the start of each row (starts with zero) 3068 . j - the column indices for each row (starts with zero) these must be sorted for each row 3069 - v - optional values in the matrix 3070 3071 Level: developer 3072 3073 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3074 3075 .keywords: matrix, aij, compressed row, sparse, sequential 3076 3077 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3078 @*/ 3079 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3080 { 3081 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 3082 PetscErrorCode ierr; 3083 3084 PetscFunctionBegin; 3085 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3086 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 3087 if (f) { 3088 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 3089 } 3090 PetscFunctionReturn(0); 3091 } 3092 3093 EXTERN_C_BEGIN 3094 #undef __FUNCT__ 3095 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3096 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3097 { 3098 PetscInt i; 3099 PetscInt m,n; 3100 PetscInt nz; 3101 PetscInt *nnz, nz_max = 0; 3102 PetscScalar *values; 3103 PetscErrorCode ierr; 3104 3105 PetscFunctionBegin; 3106 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3107 3108 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3109 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3110 for(i = 0; i < m; i++) { 3111 nz = Ii[i+1]- Ii[i]; 3112 nz_max = PetscMax(nz_max, nz); 3113 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3114 nnz[i] = nz; 3115 } 3116 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3117 ierr = PetscFree(nnz);CHKERRQ(ierr); 3118 3119 if (v) { 3120 values = (PetscScalar*) v; 3121 } else { 3122 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3123 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3124 } 3125 3126 for(i = 0; i < m; i++) { 3127 nz = Ii[i+1] - Ii[i]; 3128 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3129 } 3130 3131 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3132 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3133 3134 if (!v) { 3135 ierr = PetscFree(values);CHKERRQ(ierr); 3136 } 3137 PetscFunctionReturn(0); 3138 } 3139 EXTERN_C_END 3140 3141 #include "../src/mat/impls/dense/seq/dense.h" 3142 #include "private/petscaxpy.h" 3143 3144 #undef __FUNCT__ 3145 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3146 /* 3147 Computes (B'*A')' since computing B*A directly is untenable 3148 3149 n p p 3150 ( ) ( ) ( ) 3151 m ( A ) * n ( B ) = m ( C ) 3152 ( ) ( ) ( ) 3153 3154 */ 3155 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3156 { 3157 PetscErrorCode ierr; 3158 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3159 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3160 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3161 PetscInt i,n,m,q,p; 3162 const PetscInt *ii,*idx; 3163 const PetscScalar *b,*a,*a_q; 3164 PetscScalar *c,*c_q; 3165 3166 PetscFunctionBegin; 3167 m = A->rmap->n; 3168 n = A->cmap->n; 3169 p = B->cmap->n; 3170 a = sub_a->v; 3171 b = sub_b->a; 3172 c = sub_c->v; 3173 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3174 3175 ii = sub_b->i; 3176 idx = sub_b->j; 3177 for (i=0; i<n; i++) { 3178 q = ii[i+1] - ii[i]; 3179 while (q-->0) { 3180 c_q = c + m*(*idx); 3181 a_q = a + m*i; 3182 PetscAXPY(c_q,*b,a_q,m); 3183 idx++; 3184 b++; 3185 } 3186 } 3187 PetscFunctionReturn(0); 3188 } 3189 3190 #undef __FUNCT__ 3191 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3192 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3193 { 3194 PetscErrorCode ierr; 3195 PetscInt m=A->rmap->n,n=B->cmap->n; 3196 Mat Cmat; 3197 3198 PetscFunctionBegin; 3199 if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n); 3200 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3201 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3202 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3203 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3204 Cmat->assembled = PETSC_TRUE; 3205 *C = Cmat; 3206 PetscFunctionReturn(0); 3207 } 3208 3209 /* ----------------------------------------------------------------*/ 3210 #undef __FUNCT__ 3211 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3212 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3213 { 3214 PetscErrorCode ierr; 3215 3216 PetscFunctionBegin; 3217 if (scall == MAT_INITIAL_MATRIX){ 3218 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3219 } 3220 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3221 PetscFunctionReturn(0); 3222 } 3223 3224 3225 /*MC 3226 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3227 based on compressed sparse row format. 3228 3229 Options Database Keys: 3230 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3231 3232 Level: beginner 3233 3234 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3235 M*/ 3236 3237 EXTERN_C_BEGIN 3238 #if defined(PETSC_HAVE_PASTIX) 3239 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3240 #endif 3241 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3242 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3243 #endif 3244 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 3245 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3246 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3247 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscTruth *); 3248 #if defined(PETSC_HAVE_MUMPS) 3249 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3250 #endif 3251 #if defined(PETSC_HAVE_SUPERLU) 3252 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3253 #endif 3254 #if defined(PETSC_HAVE_SUPERLU_DIST) 3255 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3256 #endif 3257 #if defined(PETSC_HAVE_SPOOLES) 3258 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3259 #endif 3260 #if defined(PETSC_HAVE_UMFPACK) 3261 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3262 #endif 3263 #if defined(PETSC_HAVE_CHOLMOD) 3264 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3265 #endif 3266 #if defined(PETSC_HAVE_LUSOL) 3267 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3268 #endif 3269 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3270 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3271 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEnginePut_SeqAIJ(PetscObject,void*); 3272 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEngineGet_SeqAIJ(PetscObject,void*); 3273 #endif 3274 EXTERN_C_END 3275 3276 3277 EXTERN_C_BEGIN 3278 #undef __FUNCT__ 3279 #define __FUNCT__ "MatCreate_SeqAIJ" 3280 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3281 { 3282 Mat_SeqAIJ *b; 3283 PetscErrorCode ierr; 3284 PetscMPIInt size; 3285 3286 PetscFunctionBegin; 3287 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3288 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3289 3290 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3291 B->data = (void*)b; 3292 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3293 B->mapping = 0; 3294 b->row = 0; 3295 b->col = 0; 3296 b->icol = 0; 3297 b->reallocs = 0; 3298 b->ignorezeroentries = PETSC_FALSE; 3299 b->roworiented = PETSC_TRUE; 3300 b->nonew = 0; 3301 b->diag = 0; 3302 b->solve_work = 0; 3303 B->spptr = 0; 3304 b->saved_values = 0; 3305 b->idiag = 0; 3306 b->mdiag = 0; 3307 b->ssor_work = 0; 3308 b->omega = 1.0; 3309 b->fshift = 0.0; 3310 b->idiagvalid = PETSC_FALSE; 3311 b->keepnonzeropattern = PETSC_FALSE; 3312 b->xtoy = 0; 3313 b->XtoY = 0; 3314 b->compressedrow.use = PETSC_FALSE; 3315 b->compressedrow.nrows = B->rmap->n; 3316 b->compressedrow.i = PETSC_NULL; 3317 b->compressedrow.rindex = PETSC_NULL; 3318 b->compressedrow.checked = PETSC_FALSE; 3319 B->same_nonzero = PETSC_FALSE; 3320 3321 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3322 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3323 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C", 3324 "MatGetFactor_seqaij_matlab", 3325 MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3326 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3327 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3328 #endif 3329 #if defined(PETSC_HAVE_PASTIX) 3330 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 3331 "MatGetFactor_seqaij_pastix", 3332 MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3333 #endif 3334 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3335 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C", 3336 "MatGetFactor_seqaij_essl", 3337 MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3338 #endif 3339 #if defined(PETSC_HAVE_SUPERLU) 3340 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C", 3341 "MatGetFactor_seqaij_superlu", 3342 MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3343 #endif 3344 #if defined(PETSC_HAVE_SUPERLU_DIST) 3345 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C", 3346 "MatGetFactor_seqaij_superlu_dist", 3347 MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3348 #endif 3349 #if defined(PETSC_HAVE_SPOOLES) 3350 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 3351 "MatGetFactor_seqaij_spooles", 3352 MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3353 #endif 3354 #if defined(PETSC_HAVE_MUMPS) 3355 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 3356 "MatGetFactor_aij_mumps", 3357 MatGetFactor_aij_mumps);CHKERRQ(ierr); 3358 #endif 3359 #if defined(PETSC_HAVE_UMFPACK) 3360 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C", 3361 "MatGetFactor_seqaij_umfpack", 3362 MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3363 #endif 3364 #if defined(PETSC_HAVE_CHOLMOD) 3365 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 3366 "MatGetFactor_seqaij_cholmod", 3367 MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3368 #endif 3369 #if defined(PETSC_HAVE_LUSOL) 3370 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C", 3371 "MatGetFactor_seqaij_lusol", 3372 MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3373 #endif 3374 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 3375 "MatGetFactor_seqaij_petsc", 3376 MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3377 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 3378 "MatGetFactorAvailable_seqaij_petsc", 3379 MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3380 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C", 3381 "MatGetFactor_seqaij_bas", 3382 MatGetFactor_seqaij_bas); 3383 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3384 "MatSeqAIJSetColumnIndices_SeqAIJ", 3385 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3386 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3387 "MatStoreValues_SeqAIJ", 3388 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3389 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3390 "MatRetrieveValues_SeqAIJ", 3391 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3392 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3393 "MatConvert_SeqAIJ_SeqSBAIJ", 3394 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3395 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3396 "MatConvert_SeqAIJ_SeqBAIJ", 3397 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3398 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3399 "MatConvert_SeqAIJ_SeqCSRPERM", 3400 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3401 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3402 "MatConvert_SeqAIJ_SeqCRL", 3403 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3404 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3405 "MatIsTranspose_SeqAIJ", 3406 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3407 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3408 "MatIsHermitianTranspose_SeqAIJ", 3409 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3410 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3411 "MatSeqAIJSetPreallocation_SeqAIJ", 3412 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3413 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3414 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3415 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3416 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3417 "MatReorderForNonzeroDiagonal_SeqAIJ", 3418 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3419 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C", 3420 "MatMatMult_SeqDense_SeqAIJ", 3421 MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3422 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C", 3423 "MatMatMultSymbolic_SeqDense_SeqAIJ", 3424 MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3425 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C", 3426 "MatMatMultNumeric_SeqDense_SeqAIJ", 3427 MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3428 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3429 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3430 PetscFunctionReturn(0); 3431 } 3432 EXTERN_C_END 3433 3434 #undef __FUNCT__ 3435 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3436 /* 3437 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3438 */ 3439 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscTruth mallocmatspace) 3440 { 3441 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3442 PetscErrorCode ierr; 3443 PetscInt i,m = A->rmap->n; 3444 3445 PetscFunctionBegin; 3446 c = (Mat_SeqAIJ*)C->data; 3447 3448 C->factortype = A->factortype; 3449 c->row = 0; 3450 c->col = 0; 3451 c->icol = 0; 3452 c->reallocs = 0; 3453 3454 C->assembled = PETSC_TRUE; 3455 3456 ierr = PetscLayoutSetBlockSize(C->rmap,1);CHKERRQ(ierr); 3457 ierr = PetscLayoutSetBlockSize(C->cmap,1);CHKERRQ(ierr); 3458 ierr = PetscLayoutSetUp(C->rmap);CHKERRQ(ierr); 3459 ierr = PetscLayoutSetUp(C->cmap);CHKERRQ(ierr); 3460 3461 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3462 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3463 for (i=0; i<m; i++) { 3464 c->imax[i] = a->imax[i]; 3465 c->ilen[i] = a->ilen[i]; 3466 } 3467 3468 /* allocate the matrix space */ 3469 if (mallocmatspace){ 3470 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3471 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3472 c->singlemalloc = PETSC_TRUE; 3473 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3474 if (m > 0) { 3475 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3476 if (cpvalues == MAT_COPY_VALUES) { 3477 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3478 } else { 3479 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3480 } 3481 } 3482 } 3483 3484 c->ignorezeroentries = a->ignorezeroentries; 3485 c->roworiented = a->roworiented; 3486 c->nonew = a->nonew; 3487 if (a->diag) { 3488 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3489 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3490 for (i=0; i<m; i++) { 3491 c->diag[i] = a->diag[i]; 3492 } 3493 } else c->diag = 0; 3494 c->solve_work = 0; 3495 c->saved_values = 0; 3496 c->idiag = 0; 3497 c->ssor_work = 0; 3498 c->keepnonzeropattern = a->keepnonzeropattern; 3499 c->free_a = PETSC_TRUE; 3500 c->free_ij = PETSC_TRUE; 3501 c->xtoy = 0; 3502 c->XtoY = 0; 3503 3504 c->nz = a->nz; 3505 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3506 C->preallocated = PETSC_TRUE; 3507 3508 c->compressedrow.use = a->compressedrow.use; 3509 c->compressedrow.nrows = a->compressedrow.nrows; 3510 c->compressedrow.checked = a->compressedrow.checked; 3511 if (a->compressedrow.checked && a->compressedrow.use){ 3512 i = a->compressedrow.nrows; 3513 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3514 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3515 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3516 } else { 3517 c->compressedrow.use = PETSC_FALSE; 3518 c->compressedrow.i = PETSC_NULL; 3519 c->compressedrow.rindex = PETSC_NULL; 3520 } 3521 C->same_nonzero = A->same_nonzero; 3522 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3523 3524 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3525 PetscFunctionReturn(0); 3526 } 3527 3528 #undef __FUNCT__ 3529 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3530 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3531 { 3532 PetscErrorCode ierr; 3533 3534 PetscFunctionBegin; 3535 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3536 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3537 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3538 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3539 PetscFunctionReturn(0); 3540 } 3541 3542 #undef __FUNCT__ 3543 #define __FUNCT__ "MatLoad_SeqAIJ" 3544 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, const MatType type,Mat *A) 3545 { 3546 Mat_SeqAIJ *a; 3547 Mat B; 3548 PetscErrorCode ierr; 3549 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 3550 int fd; 3551 PetscMPIInt size; 3552 MPI_Comm comm; 3553 3554 PetscFunctionBegin; 3555 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3556 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3557 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3558 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3559 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3560 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3561 M = header[1]; N = header[2]; nz = header[3]; 3562 3563 if (nz < 0) { 3564 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3565 } 3566 3567 /* read in row lengths */ 3568 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3569 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3570 3571 /* check if sum of rowlengths is same as nz */ 3572 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3573 if (sum != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 3574 3575 /* create our matrix */ 3576 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 3577 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3578 ierr = MatSetType(B,type);CHKERRQ(ierr); 3579 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 3580 a = (Mat_SeqAIJ*)B->data; 3581 3582 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3583 3584 /* read in nonzero values */ 3585 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3586 3587 /* set matrix "i" values */ 3588 a->i[0] = 0; 3589 for (i=1; i<= M; i++) { 3590 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3591 a->ilen[i-1] = rowlengths[i-1]; 3592 } 3593 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3594 3595 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3596 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3597 *A = B; 3598 PetscFunctionReturn(0); 3599 } 3600 3601 #undef __FUNCT__ 3602 #define __FUNCT__ "MatLoadnew_SeqAIJ" 3603 PetscErrorCode MatLoadnew_SeqAIJ(PetscViewer viewer,Mat newMat) 3604 { 3605 Mat_SeqAIJ *a; 3606 PetscErrorCode ierr; 3607 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 3608 int fd; 3609 PetscMPIInt size; 3610 MPI_Comm comm; 3611 3612 PetscFunctionBegin; 3613 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3614 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3615 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3616 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3617 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3618 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3619 M = header[1]; N = header[2]; nz = header[3]; 3620 3621 if (nz < 0) { 3622 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3623 } 3624 3625 /* read in row lengths */ 3626 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3627 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3628 3629 /* check if sum of rowlengths is same as nz */ 3630 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3631 if (sum != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 3632 3633 /* set global size if not set already*/ 3634 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 3635 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3636 } 3637 3638 /* if sizes and type are already set, check if the vector global sizes are correct */ 3639 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 3640 if (M != rows || N != cols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols); 3641 3642 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 3643 a = (Mat_SeqAIJ*)newMat->data; 3644 3645 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3646 3647 /* read in nonzero values */ 3648 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3649 3650 /* set matrix "i" values */ 3651 a->i[0] = 0; 3652 for (i=1; i<= M; i++) { 3653 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3654 a->ilen[i-1] = rowlengths[i-1]; 3655 } 3656 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3657 3658 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3659 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3660 PetscFunctionReturn(0); 3661 } 3662 3663 #undef __FUNCT__ 3664 #define __FUNCT__ "MatEqual_SeqAIJ" 3665 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3666 { 3667 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3668 PetscErrorCode ierr; 3669 #if defined(PETSC_USE_COMPLEX) 3670 PetscInt k; 3671 #endif 3672 3673 PetscFunctionBegin; 3674 /* If the matrix dimensions are not equal,or no of nonzeros */ 3675 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 3676 *flg = PETSC_FALSE; 3677 PetscFunctionReturn(0); 3678 } 3679 3680 /* if the a->i are the same */ 3681 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3682 if (!*flg) PetscFunctionReturn(0); 3683 3684 /* if a->j are the same */ 3685 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3686 if (!*flg) PetscFunctionReturn(0); 3687 3688 /* if a->a are the same */ 3689 #if defined(PETSC_USE_COMPLEX) 3690 for (k=0; k<a->nz; k++){ 3691 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 3692 *flg = PETSC_FALSE; 3693 PetscFunctionReturn(0); 3694 } 3695 } 3696 #else 3697 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3698 #endif 3699 PetscFunctionReturn(0); 3700 } 3701 3702 #undef __FUNCT__ 3703 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3704 /*@ 3705 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3706 provided by the user. 3707 3708 Collective on MPI_Comm 3709 3710 Input Parameters: 3711 + comm - must be an MPI communicator of size 1 3712 . m - number of rows 3713 . n - number of columns 3714 . i - row indices 3715 . j - column indices 3716 - a - matrix values 3717 3718 Output Parameter: 3719 . mat - the matrix 3720 3721 Level: intermediate 3722 3723 Notes: 3724 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3725 once the matrix is destroyed 3726 3727 You cannot set new nonzero locations into this matrix, that will generate an error. 3728 3729 The i and j indices are 0 based 3730 3731 The format which is used for the sparse matrix input, is equivalent to a 3732 row-major ordering.. i.e for the following matrix, the input data expected is 3733 as shown: 3734 3735 1 0 0 3736 2 0 3 3737 4 5 6 3738 3739 i = {0,1,3,6} [size = nrow+1 = 3+1] 3740 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3741 v = {1,2,3,4,5,6} [size = nz = 6] 3742 3743 3744 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3745 3746 @*/ 3747 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3748 { 3749 PetscErrorCode ierr; 3750 PetscInt ii; 3751 Mat_SeqAIJ *aij; 3752 #if defined(PETSC_USE_DEBUG) 3753 PetscInt jj; 3754 #endif 3755 3756 PetscFunctionBegin; 3757 if (i[0]) { 3758 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3759 } 3760 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3761 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3762 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3763 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3764 aij = (Mat_SeqAIJ*)(*mat)->data; 3765 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3766 3767 aij->i = i; 3768 aij->j = j; 3769 aij->a = a; 3770 aij->singlemalloc = PETSC_FALSE; 3771 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3772 aij->free_a = PETSC_FALSE; 3773 aij->free_ij = PETSC_FALSE; 3774 3775 for (ii=0; ii<m; ii++) { 3776 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3777 #if defined(PETSC_USE_DEBUG) 3778 if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 3779 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3780 if (j[jj] < j[jj-1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is not sorted",jj-i[ii],j[jj],ii); 3781 if (j[jj] == j[jj]-1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is identical to previous entry",jj-i[ii],j[jj],ii); 3782 } 3783 #endif 3784 } 3785 #if defined(PETSC_USE_DEBUG) 3786 for (ii=0; ii<aij->i[m]; ii++) { 3787 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3788 if (j[ii] > n - 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3789 } 3790 #endif 3791 3792 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3793 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3794 PetscFunctionReturn(0); 3795 } 3796 3797 #undef __FUNCT__ 3798 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3799 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3800 { 3801 PetscErrorCode ierr; 3802 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3803 3804 PetscFunctionBegin; 3805 if (coloring->ctype == IS_COLORING_GLOBAL) { 3806 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3807 a->coloring = coloring; 3808 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3809 PetscInt i,*larray; 3810 ISColoring ocoloring; 3811 ISColoringValue *colors; 3812 3813 /* set coloring for diagonal portion */ 3814 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3815 for (i=0; i<A->cmap->n; i++) { 3816 larray[i] = i; 3817 } 3818 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3819 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3820 for (i=0; i<A->cmap->n; i++) { 3821 colors[i] = coloring->colors[larray[i]]; 3822 } 3823 ierr = PetscFree(larray);CHKERRQ(ierr); 3824 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 3825 a->coloring = ocoloring; 3826 } 3827 PetscFunctionReturn(0); 3828 } 3829 3830 #if defined(PETSC_HAVE_ADIC) 3831 EXTERN_C_BEGIN 3832 #include "adic/ad_utils.h" 3833 EXTERN_C_END 3834 3835 #undef __FUNCT__ 3836 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3837 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3838 { 3839 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3840 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3841 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3842 ISColoringValue *color; 3843 3844 PetscFunctionBegin; 3845 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3846 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3847 color = a->coloring->colors; 3848 /* loop over rows */ 3849 for (i=0; i<m; i++) { 3850 nz = ii[i+1] - ii[i]; 3851 /* loop over columns putting computed value into matrix */ 3852 for (j=0; j<nz; j++) { 3853 *v++ = values[color[*jj++]]; 3854 } 3855 values += nlen; /* jump to next row of derivatives */ 3856 } 3857 PetscFunctionReturn(0); 3858 } 3859 #endif 3860 3861 #undef __FUNCT__ 3862 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3863 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3864 { 3865 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3866 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 3867 MatScalar *v = a->a; 3868 PetscScalar *values = (PetscScalar *)advalues; 3869 ISColoringValue *color; 3870 3871 PetscFunctionBegin; 3872 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3873 color = a->coloring->colors; 3874 /* loop over rows */ 3875 for (i=0; i<m; i++) { 3876 nz = ii[i+1] - ii[i]; 3877 /* loop over columns putting computed value into matrix */ 3878 for (j=0; j<nz; j++) { 3879 *v++ = values[color[*jj++]]; 3880 } 3881 values += nl; /* jump to next row of derivatives */ 3882 } 3883 PetscFunctionReturn(0); 3884 } 3885 3886 /* 3887 Special version for direct calls from Fortran 3888 */ 3889 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3890 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3891 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3892 #define matsetvaluesseqaij_ matsetvaluesseqaij 3893 #endif 3894 3895 /* Change these macros so can be used in void function */ 3896 #undef CHKERRQ 3897 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3898 #undef SETERRQ2 3899 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 3900 3901 EXTERN_C_BEGIN 3902 #undef __FUNCT__ 3903 #define __FUNCT__ "matsetvaluesseqaij_" 3904 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3905 { 3906 Mat A = *AA; 3907 PetscInt m = *mm, n = *nn; 3908 InsertMode is = *isis; 3909 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3910 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3911 PetscInt *imax,*ai,*ailen; 3912 PetscErrorCode ierr; 3913 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3914 MatScalar *ap,value,*aa; 3915 PetscTruth ignorezeroentries = a->ignorezeroentries; 3916 PetscTruth roworiented = a->roworiented; 3917 3918 PetscFunctionBegin; 3919 ierr = MatPreallocated(A);CHKERRQ(ierr); 3920 imax = a->imax; 3921 ai = a->i; 3922 ailen = a->ilen; 3923 aj = a->j; 3924 aa = a->a; 3925 3926 for (k=0; k<m; k++) { /* loop over added rows */ 3927 row = im[k]; 3928 if (row < 0) continue; 3929 #if defined(PETSC_USE_DEBUG) 3930 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3931 #endif 3932 rp = aj + ai[row]; ap = aa + ai[row]; 3933 rmax = imax[row]; nrow = ailen[row]; 3934 low = 0; 3935 high = nrow; 3936 for (l=0; l<n; l++) { /* loop over added columns */ 3937 if (in[l] < 0) continue; 3938 #if defined(PETSC_USE_DEBUG) 3939 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3940 #endif 3941 col = in[l]; 3942 if (roworiented) { 3943 value = v[l + k*n]; 3944 } else { 3945 value = v[k + l*m]; 3946 } 3947 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3948 3949 if (col <= lastcol) low = 0; else high = nrow; 3950 lastcol = col; 3951 while (high-low > 5) { 3952 t = (low+high)/2; 3953 if (rp[t] > col) high = t; 3954 else low = t; 3955 } 3956 for (i=low; i<high; i++) { 3957 if (rp[i] > col) break; 3958 if (rp[i] == col) { 3959 if (is == ADD_VALUES) ap[i] += value; 3960 else ap[i] = value; 3961 goto noinsert; 3962 } 3963 } 3964 if (value == 0.0 && ignorezeroentries) goto noinsert; 3965 if (nonew == 1) goto noinsert; 3966 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3967 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3968 N = nrow++ - 1; a->nz++; high++; 3969 /* shift up all the later entries in this row */ 3970 for (ii=N; ii>=i; ii--) { 3971 rp[ii+1] = rp[ii]; 3972 ap[ii+1] = ap[ii]; 3973 } 3974 rp[i] = col; 3975 ap[i] = value; 3976 noinsert:; 3977 low = i + 1; 3978 } 3979 ailen[row] = nrow; 3980 } 3981 A->same_nonzero = PETSC_FALSE; 3982 PetscFunctionReturnVoid(); 3983 } 3984 EXTERN_C_END 3985