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 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1004 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1005 aj = a->j; 1006 aa = a->a; 1007 ii = a->i; 1008 if (usecprow){ /* use compressed row format */ 1009 m = a->compressedrow.nrows; 1010 ii = a->compressedrow.i; 1011 ridx = a->compressedrow.rindex; 1012 for (i=0; i<m; i++){ 1013 n = ii[i+1] - ii[i]; 1014 aj = a->j + ii[i]; 1015 aa = a->a + ii[i]; 1016 sum = 0.0; 1017 nonzerorow += (n>0); 1018 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1019 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1020 y[*ridx++] = sum; 1021 } 1022 } else { /* do not use compressed row format */ 1023 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1024 fortranmultaij_(&m,x,ii,aj,aa,y); 1025 #else 1026 #if defined(PETSC_HAVE_CUDA) 1027 ierr = VecCUDACopyToGPU(xx);CHKERRQ(ierr); 1028 cusp::multiply(a->GPUmatrix,xx->GPUarray,yy->GPUarray); 1029 ierr = PetscPrintf(PETSC_COMM_WORLD,"MULTIPLIED WITH CUSP");CHKERRQ(ierr); 1030 yy->valid_GPU_array = PETSC_CUDA_GPU; 1031 #else 1032 for (i=0; i<m; i++) { 1033 n = ii[i+1] - ii[i]; 1034 aj = a->j + ii[i]; 1035 aa = a->a + ii[i]; 1036 sum = 0.0; 1037 nonzerorow += (n>0); 1038 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1039 y[i] = sum; 1040 } 1041 #endif 1042 #endif 1043 } 1044 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1045 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1046 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1047 PetscFunctionReturn(0); 1048 } 1049 1050 #include "../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h" 1051 #undef __FUNCT__ 1052 #define __FUNCT__ "MatMultAdd_SeqAIJ" 1053 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1054 { 1055 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1056 PetscScalar *x,*y,*z; 1057 const MatScalar *aa; 1058 PetscErrorCode ierr; 1059 PetscInt m = A->rmap->n,*aj,*ii; 1060 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1061 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 1062 PetscScalar sum; 1063 PetscTruth usecprow=a->compressedrow.use; 1064 #endif 1065 1066 PetscFunctionBegin; 1067 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1068 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1069 if (zz != yy) { 1070 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1071 } else { 1072 z = y; 1073 } 1074 1075 aj = a->j; 1076 aa = a->a; 1077 ii = a->i; 1078 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1079 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1080 #else 1081 if (usecprow){ /* use compressed row format */ 1082 if (zz != yy){ 1083 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1084 } 1085 m = a->compressedrow.nrows; 1086 ii = a->compressedrow.i; 1087 ridx = a->compressedrow.rindex; 1088 for (i=0; i<m; i++){ 1089 n = ii[i+1] - ii[i]; 1090 aj = a->j + ii[i]; 1091 aa = a->a + ii[i]; 1092 sum = y[*ridx]; 1093 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 1094 z[*ridx++] = sum; 1095 } 1096 } else { /* do not use compressed row format */ 1097 for (i=0; i<m; i++) { 1098 jrow = ii[i]; 1099 n = ii[i+1] - jrow; 1100 sum = y[i]; 1101 for (j=0; j<n; j++) { 1102 sum += aa[jrow]*x[aj[jrow]]; jrow++; 1103 } 1104 z[i] = sum; 1105 } 1106 } 1107 #endif 1108 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1109 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1110 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1111 if (zz != yy) { 1112 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1113 } 1114 PetscFunctionReturn(0); 1115 } 1116 1117 /* 1118 Adds diagonal pointers to sparse matrix structure. 1119 */ 1120 #undef __FUNCT__ 1121 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1122 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1123 { 1124 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1125 PetscErrorCode ierr; 1126 PetscInt i,j,m = A->rmap->n; 1127 1128 PetscFunctionBegin; 1129 if (!a->diag) { 1130 ierr = PetscMalloc(m*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 1131 ierr = PetscLogObjectMemory(A, m*sizeof(PetscInt));CHKERRQ(ierr); 1132 } 1133 for (i=0; i<A->rmap->n; i++) { 1134 a->diag[i] = a->i[i+1]; 1135 for (j=a->i[i]; j<a->i[i+1]; j++) { 1136 if (a->j[j] == i) { 1137 a->diag[i] = j; 1138 break; 1139 } 1140 } 1141 } 1142 PetscFunctionReturn(0); 1143 } 1144 1145 /* 1146 Checks for missing diagonals 1147 */ 1148 #undef __FUNCT__ 1149 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1150 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscTruth *missing,PetscInt *d) 1151 { 1152 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1153 PetscInt *diag,*jj = a->j,i; 1154 1155 PetscFunctionBegin; 1156 *missing = PETSC_FALSE; 1157 if (A->rmap->n > 0 && !jj) { 1158 *missing = PETSC_TRUE; 1159 if (d) *d = 0; 1160 PetscInfo(A,"Matrix has no entries therefor is missing diagonal"); 1161 } else { 1162 diag = a->diag; 1163 for (i=0; i<A->rmap->n; i++) { 1164 if (jj[diag[i]] != i) { 1165 *missing = PETSC_TRUE; 1166 if (d) *d = i; 1167 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1168 } 1169 } 1170 } 1171 PetscFunctionReturn(0); 1172 } 1173 1174 EXTERN_C_BEGIN 1175 #undef __FUNCT__ 1176 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1177 PetscErrorCode PETSCMAT_DLLEXPORT MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1178 { 1179 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1180 PetscErrorCode ierr; 1181 PetscInt i,*diag,m = A->rmap->n; 1182 MatScalar *v = a->a; 1183 PetscScalar *idiag,*mdiag; 1184 1185 PetscFunctionBegin; 1186 if (a->idiagvalid) PetscFunctionReturn(0); 1187 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1188 diag = a->diag; 1189 if (!a->idiag) { 1190 ierr = PetscMalloc3(m,PetscScalar,&a->idiag,m,PetscScalar,&a->mdiag,m,PetscScalar,&a->ssor_work);CHKERRQ(ierr); 1191 ierr = PetscLogObjectMemory(A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1192 v = a->a; 1193 } 1194 mdiag = a->mdiag; 1195 idiag = a->idiag; 1196 1197 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1198 for (i=0; i<m; i++) { 1199 mdiag[i] = v[diag[i]]; 1200 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1201 idiag[i] = 1.0/v[diag[i]]; 1202 } 1203 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1204 } else { 1205 for (i=0; i<m; i++) { 1206 mdiag[i] = v[diag[i]]; 1207 idiag[i] = omega/(fshift + v[diag[i]]); 1208 } 1209 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1210 } 1211 a->idiagvalid = PETSC_TRUE; 1212 PetscFunctionReturn(0); 1213 } 1214 EXTERN_C_END 1215 1216 #include "../src/mat/impls/aij/seq/ftn-kernels/frelax.h" 1217 #undef __FUNCT__ 1218 #define __FUNCT__ "MatSOR_SeqAIJ" 1219 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1220 { 1221 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1222 PetscScalar *x,d,sum,*t,scale; 1223 const MatScalar *v = a->a,*idiag=0,*mdiag; 1224 const PetscScalar *b, *bs,*xb, *ts; 1225 PetscErrorCode ierr; 1226 PetscInt n = A->cmap->n,m = A->rmap->n,i; 1227 const PetscInt *idx,*diag; 1228 1229 PetscFunctionBegin; 1230 its = its*lits; 1231 1232 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1233 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1234 a->fshift = fshift; 1235 a->omega = omega; 1236 1237 diag = a->diag; 1238 t = a->ssor_work; 1239 idiag = a->idiag; 1240 mdiag = a->mdiag; 1241 1242 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1243 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1244 CHKMEMQ; 1245 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1246 if (flag == SOR_APPLY_UPPER) { 1247 /* apply (U + D/omega) to the vector */ 1248 bs = b; 1249 for (i=0; i<m; i++) { 1250 d = fshift + mdiag[i]; 1251 n = a->i[i+1] - diag[i] - 1; 1252 idx = a->j + diag[i] + 1; 1253 v = a->a + diag[i] + 1; 1254 sum = b[i]*d/omega; 1255 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1256 x[i] = sum; 1257 } 1258 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1259 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1260 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1261 PetscFunctionReturn(0); 1262 } 1263 1264 if (flag == SOR_APPLY_LOWER) { 1265 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1266 } else if (flag & SOR_EISENSTAT) { 1267 /* Let A = L + U + D; where L is lower trianglar, 1268 U is upper triangular, E = D/omega; This routine applies 1269 1270 (L + E)^{-1} A (U + E)^{-1} 1271 1272 to a vector efficiently using Eisenstat's trick. 1273 */ 1274 scale = (2.0/omega) - 1.0; 1275 1276 /* x = (E + U)^{-1} b */ 1277 for (i=m-1; i>=0; i--) { 1278 n = a->i[i+1] - diag[i] - 1; 1279 idx = a->j + diag[i] + 1; 1280 v = a->a + diag[i] + 1; 1281 sum = b[i]; 1282 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1283 x[i] = sum*idiag[i]; 1284 } 1285 1286 /* t = b - (2*E - D)x */ 1287 v = a->a; 1288 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1289 1290 /* t = (E + L)^{-1}t */ 1291 ts = t; 1292 diag = a->diag; 1293 for (i=0; i<m; i++) { 1294 n = diag[i] - a->i[i]; 1295 idx = a->j + a->i[i]; 1296 v = a->a + a->i[i]; 1297 sum = t[i]; 1298 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1299 t[i] = sum*idiag[i]; 1300 /* x = x + t */ 1301 x[i] += t[i]; 1302 } 1303 1304 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1305 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1306 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1307 PetscFunctionReturn(0); 1308 } 1309 if (flag & SOR_ZERO_INITIAL_GUESS) { 1310 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1311 for (i=0; i<m; i++) { 1312 n = diag[i] - a->i[i]; 1313 idx = a->j + a->i[i]; 1314 v = a->a + a->i[i]; 1315 sum = b[i]; 1316 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1317 t[i] = sum; 1318 x[i] = sum*idiag[i]; 1319 } 1320 xb = t; 1321 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1322 } else xb = b; 1323 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1324 for (i=m-1; i>=0; i--) { 1325 n = a->i[i+1] - diag[i] - 1; 1326 idx = a->j + diag[i] + 1; 1327 v = a->a + diag[i] + 1; 1328 sum = xb[i]; 1329 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1330 if (xb == b) { 1331 x[i] = sum*idiag[i]; 1332 } else { 1333 x[i] = (1-omega)*x[i] + sum*idiag[i]; 1334 } 1335 } 1336 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1337 } 1338 its--; 1339 } 1340 while (its--) { 1341 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1342 for (i=0; i<m; i++) { 1343 n = a->i[i+1] - a->i[i]; 1344 idx = a->j + a->i[i]; 1345 v = a->a + a->i[i]; 1346 sum = b[i]; 1347 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1348 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1349 } 1350 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1351 } 1352 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1353 for (i=m-1; i>=0; i--) { 1354 n = a->i[i+1] - a->i[i]; 1355 idx = a->j + a->i[i]; 1356 v = a->a + a->i[i]; 1357 sum = b[i]; 1358 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1359 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1360 } 1361 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1362 } 1363 } 1364 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1365 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1366 CHKMEMQ; PetscFunctionReturn(0); 1367 } 1368 1369 1370 #undef __FUNCT__ 1371 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1372 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1373 { 1374 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1375 1376 PetscFunctionBegin; 1377 info->block_size = 1.0; 1378 info->nz_allocated = (double)a->maxnz; 1379 info->nz_used = (double)a->nz; 1380 info->nz_unneeded = (double)(a->maxnz - a->nz); 1381 info->assemblies = (double)A->num_ass; 1382 info->mallocs = (double)A->info.mallocs; 1383 info->memory = ((PetscObject)A)->mem; 1384 if (A->factortype) { 1385 info->fill_ratio_given = A->info.fill_ratio_given; 1386 info->fill_ratio_needed = A->info.fill_ratio_needed; 1387 info->factor_mallocs = A->info.factor_mallocs; 1388 } else { 1389 info->fill_ratio_given = 0; 1390 info->fill_ratio_needed = 0; 1391 info->factor_mallocs = 0; 1392 } 1393 PetscFunctionReturn(0); 1394 } 1395 1396 #undef __FUNCT__ 1397 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1398 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1399 { 1400 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1401 PetscInt i,m = A->rmap->n - 1,d = 0; 1402 PetscErrorCode ierr; 1403 PetscTruth missing; 1404 1405 PetscFunctionBegin; 1406 if (a->keepnonzeropattern) { 1407 for (i=0; i<N; i++) { 1408 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1409 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1410 } 1411 if (diag != 0.0) { 1412 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1413 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1414 for (i=0; i<N; i++) { 1415 a->a[a->diag[rows[i]]] = diag; 1416 } 1417 } 1418 A->same_nonzero = PETSC_TRUE; 1419 } else { 1420 if (diag != 0.0) { 1421 for (i=0; i<N; i++) { 1422 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1423 if (a->ilen[rows[i]] > 0) { 1424 a->ilen[rows[i]] = 1; 1425 a->a[a->i[rows[i]]] = diag; 1426 a->j[a->i[rows[i]]] = rows[i]; 1427 } else { /* in case row was completely empty */ 1428 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1429 } 1430 } 1431 } else { 1432 for (i=0; i<N; i++) { 1433 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1434 a->ilen[rows[i]] = 0; 1435 } 1436 } 1437 A->same_nonzero = PETSC_FALSE; 1438 } 1439 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1440 PetscFunctionReturn(0); 1441 } 1442 1443 #undef __FUNCT__ 1444 #define __FUNCT__ "MatGetRow_SeqAIJ" 1445 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1446 { 1447 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1448 PetscInt *itmp; 1449 1450 PetscFunctionBegin; 1451 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1452 1453 *nz = a->i[row+1] - a->i[row]; 1454 if (v) *v = a->a + a->i[row]; 1455 if (idx) { 1456 itmp = a->j + a->i[row]; 1457 if (*nz) { 1458 *idx = itmp; 1459 } 1460 else *idx = 0; 1461 } 1462 PetscFunctionReturn(0); 1463 } 1464 1465 /* remove this function? */ 1466 #undef __FUNCT__ 1467 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1468 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1469 { 1470 PetscFunctionBegin; 1471 PetscFunctionReturn(0); 1472 } 1473 1474 #undef __FUNCT__ 1475 #define __FUNCT__ "MatNorm_SeqAIJ" 1476 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1477 { 1478 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1479 MatScalar *v = a->a; 1480 PetscReal sum = 0.0; 1481 PetscErrorCode ierr; 1482 PetscInt i,j; 1483 1484 PetscFunctionBegin; 1485 if (type == NORM_FROBENIUS) { 1486 for (i=0; i<a->nz; i++) { 1487 #if defined(PETSC_USE_COMPLEX) 1488 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1489 #else 1490 sum += (*v)*(*v); v++; 1491 #endif 1492 } 1493 *nrm = sqrt(sum); 1494 } else if (type == NORM_1) { 1495 PetscReal *tmp; 1496 PetscInt *jj = a->j; 1497 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1498 ierr = PetscMemzero(tmp,A->cmap->n*sizeof(PetscReal));CHKERRQ(ierr); 1499 *nrm = 0.0; 1500 for (j=0; j<a->nz; j++) { 1501 tmp[*jj++] += PetscAbsScalar(*v); v++; 1502 } 1503 for (j=0; j<A->cmap->n; j++) { 1504 if (tmp[j] > *nrm) *nrm = tmp[j]; 1505 } 1506 ierr = PetscFree(tmp);CHKERRQ(ierr); 1507 } else if (type == NORM_INFINITY) { 1508 *nrm = 0.0; 1509 for (j=0; j<A->rmap->n; j++) { 1510 v = a->a + a->i[j]; 1511 sum = 0.0; 1512 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1513 sum += PetscAbsScalar(*v); v++; 1514 } 1515 if (sum > *nrm) *nrm = sum; 1516 } 1517 } else { 1518 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1519 } 1520 PetscFunctionReturn(0); 1521 } 1522 1523 #undef __FUNCT__ 1524 #define __FUNCT__ "MatTranspose_SeqAIJ" 1525 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 1526 { 1527 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1528 Mat C; 1529 PetscErrorCode ierr; 1530 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 1531 MatScalar *array = a->a; 1532 1533 PetscFunctionBegin; 1534 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"); 1535 1536 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 1537 ierr = PetscMalloc((1+A->cmap->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1538 ierr = PetscMemzero(col,(1+A->cmap->n)*sizeof(PetscInt));CHKERRQ(ierr); 1539 1540 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1541 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1542 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 1543 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1544 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1545 ierr = PetscFree(col);CHKERRQ(ierr); 1546 } else { 1547 C = *B; 1548 } 1549 1550 for (i=0; i<m; i++) { 1551 len = ai[i+1]-ai[i]; 1552 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1553 array += len; 1554 aj += len; 1555 } 1556 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1557 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1558 1559 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 1560 *B = C; 1561 } else { 1562 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1563 } 1564 PetscFunctionReturn(0); 1565 } 1566 1567 EXTERN_C_BEGIN 1568 #undef __FUNCT__ 1569 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1570 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1571 { 1572 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1573 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1574 MatScalar *va,*vb; 1575 PetscErrorCode ierr; 1576 PetscInt ma,na,mb,nb, i; 1577 1578 PetscFunctionBegin; 1579 bij = (Mat_SeqAIJ *) B->data; 1580 1581 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1582 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1583 if (ma!=nb || na!=mb){ 1584 *f = PETSC_FALSE; 1585 PetscFunctionReturn(0); 1586 } 1587 aii = aij->i; bii = bij->i; 1588 adx = aij->j; bdx = bij->j; 1589 va = aij->a; vb = bij->a; 1590 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1591 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1592 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1593 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1594 1595 *f = PETSC_TRUE; 1596 for (i=0; i<ma; i++) { 1597 while (aptr[i]<aii[i+1]) { 1598 PetscInt idc,idr; 1599 PetscScalar vc,vr; 1600 /* column/row index/value */ 1601 idc = adx[aptr[i]]; 1602 idr = bdx[bptr[idc]]; 1603 vc = va[aptr[i]]; 1604 vr = vb[bptr[idc]]; 1605 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1606 *f = PETSC_FALSE; 1607 goto done; 1608 } else { 1609 aptr[i]++; 1610 if (B || i!=idc) bptr[idc]++; 1611 } 1612 } 1613 } 1614 done: 1615 ierr = PetscFree(aptr);CHKERRQ(ierr); 1616 if (B) { 1617 ierr = PetscFree(bptr);CHKERRQ(ierr); 1618 } 1619 PetscFunctionReturn(0); 1620 } 1621 EXTERN_C_END 1622 1623 EXTERN_C_BEGIN 1624 #undef __FUNCT__ 1625 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 1626 PetscErrorCode PETSCMAT_DLLEXPORT MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1627 { 1628 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1629 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1630 MatScalar *va,*vb; 1631 PetscErrorCode ierr; 1632 PetscInt ma,na,mb,nb, i; 1633 1634 PetscFunctionBegin; 1635 bij = (Mat_SeqAIJ *) B->data; 1636 1637 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1638 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1639 if (ma!=nb || na!=mb){ 1640 *f = PETSC_FALSE; 1641 PetscFunctionReturn(0); 1642 } 1643 aii = aij->i; bii = bij->i; 1644 adx = aij->j; bdx = bij->j; 1645 va = aij->a; vb = bij->a; 1646 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1647 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1648 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1649 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1650 1651 *f = PETSC_TRUE; 1652 for (i=0; i<ma; i++) { 1653 while (aptr[i]<aii[i+1]) { 1654 PetscInt idc,idr; 1655 PetscScalar vc,vr; 1656 /* column/row index/value */ 1657 idc = adx[aptr[i]]; 1658 idr = bdx[bptr[idc]]; 1659 vc = va[aptr[i]]; 1660 vr = vb[bptr[idc]]; 1661 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 1662 *f = PETSC_FALSE; 1663 goto done; 1664 } else { 1665 aptr[i]++; 1666 if (B || i!=idc) bptr[idc]++; 1667 } 1668 } 1669 } 1670 done: 1671 ierr = PetscFree(aptr);CHKERRQ(ierr); 1672 if (B) { 1673 ierr = PetscFree(bptr);CHKERRQ(ierr); 1674 } 1675 PetscFunctionReturn(0); 1676 } 1677 EXTERN_C_END 1678 1679 #undef __FUNCT__ 1680 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1681 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1682 { 1683 PetscErrorCode ierr; 1684 PetscFunctionBegin; 1685 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1686 PetscFunctionReturn(0); 1687 } 1688 1689 #undef __FUNCT__ 1690 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 1691 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1692 { 1693 PetscErrorCode ierr; 1694 PetscFunctionBegin; 1695 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1696 PetscFunctionReturn(0); 1697 } 1698 1699 #undef __FUNCT__ 1700 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1701 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1702 { 1703 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1704 PetscScalar *l,*r,x; 1705 MatScalar *v; 1706 PetscErrorCode ierr; 1707 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 1708 1709 PetscFunctionBegin; 1710 if (ll) { 1711 /* The local size is used so that VecMPI can be passed to this routine 1712 by MatDiagonalScale_MPIAIJ */ 1713 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1714 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1715 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1716 v = a->a; 1717 for (i=0; i<m; i++) { 1718 x = l[i]; 1719 M = a->i[i+1] - a->i[i]; 1720 for (j=0; j<M; j++) { (*v++) *= x;} 1721 } 1722 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1723 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1724 } 1725 if (rr) { 1726 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1727 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1728 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1729 v = a->a; jj = a->j; 1730 for (i=0; i<nz; i++) { 1731 (*v++) *= r[*jj++]; 1732 } 1733 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1734 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1735 } 1736 PetscFunctionReturn(0); 1737 } 1738 1739 #undef __FUNCT__ 1740 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1741 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1742 { 1743 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1744 PetscErrorCode ierr; 1745 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 1746 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1747 const PetscInt *irow,*icol; 1748 PetscInt nrows,ncols; 1749 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1750 MatScalar *a_new,*mat_a; 1751 Mat C; 1752 PetscTruth stride,sorted; 1753 1754 PetscFunctionBegin; 1755 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 1756 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1757 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 1758 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1759 1760 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1761 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1762 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1763 1764 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1765 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1766 if (stride && step == 1) { 1767 /* special case of contiguous rows */ 1768 ierr = PetscMalloc2(nrows,PetscInt,&lens,nrows,PetscInt,&starts);CHKERRQ(ierr); 1769 /* loop over new rows determining lens and starting points */ 1770 for (i=0; i<nrows; i++) { 1771 kstart = ai[irow[i]]; 1772 kend = kstart + ailen[irow[i]]; 1773 for (k=kstart; k<kend; k++) { 1774 if (aj[k] >= first) { 1775 starts[i] = k; 1776 break; 1777 } 1778 } 1779 sum = 0; 1780 while (k < kend) { 1781 if (aj[k++] >= first+ncols) break; 1782 sum++; 1783 } 1784 lens[i] = sum; 1785 } 1786 /* create submatrix */ 1787 if (scall == MAT_REUSE_MATRIX) { 1788 PetscInt n_cols,n_rows; 1789 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1790 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1791 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1792 C = *B; 1793 } else { 1794 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1795 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1796 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1797 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1798 } 1799 c = (Mat_SeqAIJ*)C->data; 1800 1801 /* loop over rows inserting into submatrix */ 1802 a_new = c->a; 1803 j_new = c->j; 1804 i_new = c->i; 1805 1806 for (i=0; i<nrows; i++) { 1807 ii = starts[i]; 1808 lensi = lens[i]; 1809 for (k=0; k<lensi; k++) { 1810 *j_new++ = aj[ii+k] - first; 1811 } 1812 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1813 a_new += lensi; 1814 i_new[i+1] = i_new[i] + lensi; 1815 c->ilen[i] = lensi; 1816 } 1817 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 1818 } else { 1819 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1820 ierr = PetscMalloc(oldcols*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1821 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1822 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1823 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1824 /* determine lens of each row */ 1825 for (i=0; i<nrows; i++) { 1826 kstart = ai[irow[i]]; 1827 kend = kstart + a->ilen[irow[i]]; 1828 lens[i] = 0; 1829 for (k=kstart; k<kend; k++) { 1830 if (smap[aj[k]]) { 1831 lens[i]++; 1832 } 1833 } 1834 } 1835 /* Create and fill new matrix */ 1836 if (scall == MAT_REUSE_MATRIX) { 1837 PetscTruth equal; 1838 1839 c = (Mat_SeqAIJ *)((*B)->data); 1840 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1841 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1842 if (!equal) { 1843 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1844 } 1845 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 1846 C = *B; 1847 } else { 1848 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1849 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1850 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1851 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1852 } 1853 c = (Mat_SeqAIJ *)(C->data); 1854 for (i=0; i<nrows; i++) { 1855 row = irow[i]; 1856 kstart = ai[row]; 1857 kend = kstart + a->ilen[row]; 1858 mat_i = c->i[i]; 1859 mat_j = c->j + mat_i; 1860 mat_a = c->a + mat_i; 1861 mat_ilen = c->ilen + i; 1862 for (k=kstart; k<kend; k++) { 1863 if ((tcol=smap[a->j[k]])) { 1864 *mat_j++ = tcol - 1; 1865 *mat_a++ = a->a[k]; 1866 (*mat_ilen)++; 1867 1868 } 1869 } 1870 } 1871 /* Free work space */ 1872 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1873 ierr = PetscFree(smap);CHKERRQ(ierr); 1874 ierr = PetscFree(lens);CHKERRQ(ierr); 1875 } 1876 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1877 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1878 1879 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1880 *B = C; 1881 PetscFunctionReturn(0); 1882 } 1883 1884 #undef __FUNCT__ 1885 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1886 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 1887 { 1888 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1889 PetscErrorCode ierr; 1890 Mat outA; 1891 PetscTruth row_identity,col_identity; 1892 1893 PetscFunctionBegin; 1894 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1895 1896 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1897 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1898 1899 outA = inA; 1900 outA->factortype = MAT_FACTOR_LU; 1901 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1902 if (a->row) { ierr = ISDestroy(a->row);CHKERRQ(ierr);} 1903 a->row = row; 1904 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1905 if (a->col) { ierr = ISDestroy(a->col);CHKERRQ(ierr);} 1906 a->col = col; 1907 1908 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1909 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1910 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1911 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1912 1913 if (!a->solve_work) { /* this matrix may have been factored before */ 1914 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1915 ierr = PetscLogObjectMemory(inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 1916 } 1917 1918 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1919 if (row_identity && col_identity) { 1920 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 1921 } else { 1922 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 1923 } 1924 PetscFunctionReturn(0); 1925 } 1926 1927 #undef __FUNCT__ 1928 #define __FUNCT__ "MatScale_SeqAIJ" 1929 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1930 { 1931 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1932 PetscScalar oalpha = alpha; 1933 PetscErrorCode ierr; 1934 PetscBLASInt one = 1,bnz = PetscBLASIntCast(a->nz); 1935 1936 PetscFunctionBegin; 1937 BLASscal_(&bnz,&oalpha,a->a,&one); 1938 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1939 PetscFunctionReturn(0); 1940 } 1941 1942 #undef __FUNCT__ 1943 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1944 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1945 { 1946 PetscErrorCode ierr; 1947 PetscInt i; 1948 1949 PetscFunctionBegin; 1950 if (scall == MAT_INITIAL_MATRIX) { 1951 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1952 } 1953 1954 for (i=0; i<n; i++) { 1955 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1956 } 1957 PetscFunctionReturn(0); 1958 } 1959 1960 #undef __FUNCT__ 1961 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1962 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1963 { 1964 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1965 PetscErrorCode ierr; 1966 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 1967 const PetscInt *idx; 1968 PetscInt start,end,*ai,*aj; 1969 PetscBT table; 1970 1971 PetscFunctionBegin; 1972 m = A->rmap->n; 1973 ai = a->i; 1974 aj = a->j; 1975 1976 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1977 1978 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1979 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1980 1981 for (i=0; i<is_max; i++) { 1982 /* Initialize the two local arrays */ 1983 isz = 0; 1984 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1985 1986 /* Extract the indices, assume there can be duplicate entries */ 1987 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1988 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1989 1990 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1991 for (j=0; j<n ; ++j){ 1992 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1993 } 1994 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1995 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1996 1997 k = 0; 1998 for (j=0; j<ov; j++){ /* for each overlap */ 1999 n = isz; 2000 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 2001 row = nidx[k]; 2002 start = ai[row]; 2003 end = ai[row+1]; 2004 for (l = start; l<end ; l++){ 2005 val = aj[l] ; 2006 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 2007 } 2008 } 2009 } 2010 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 2011 } 2012 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 2013 ierr = PetscFree(nidx);CHKERRQ(ierr); 2014 PetscFunctionReturn(0); 2015 } 2016 2017 /* -------------------------------------------------------------- */ 2018 #undef __FUNCT__ 2019 #define __FUNCT__ "MatPermute_SeqAIJ" 2020 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2021 { 2022 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2023 PetscErrorCode ierr; 2024 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2025 const PetscInt *row,*col; 2026 PetscInt *cnew,j,*lens; 2027 IS icolp,irowp; 2028 PetscInt *cwork = PETSC_NULL; 2029 PetscScalar *vwork = PETSC_NULL; 2030 2031 PetscFunctionBegin; 2032 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2033 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2034 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2035 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2036 2037 /* determine lengths of permuted rows */ 2038 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2039 for (i=0; i<m; i++) { 2040 lens[row[i]] = a->i[i+1] - a->i[i]; 2041 } 2042 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 2043 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2044 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2045 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2046 ierr = PetscFree(lens);CHKERRQ(ierr); 2047 2048 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2049 for (i=0; i<m; i++) { 2050 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2051 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 2052 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2053 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2054 } 2055 ierr = PetscFree(cnew);CHKERRQ(ierr); 2056 (*B)->assembled = PETSC_FALSE; 2057 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2058 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2059 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2060 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2061 ierr = ISDestroy(irowp);CHKERRQ(ierr); 2062 ierr = ISDestroy(icolp);CHKERRQ(ierr); 2063 PetscFunctionReturn(0); 2064 } 2065 2066 #undef __FUNCT__ 2067 #define __FUNCT__ "MatCopy_SeqAIJ" 2068 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2069 { 2070 PetscErrorCode ierr; 2071 2072 PetscFunctionBegin; 2073 /* If the two matrices have the same copy implementation, use fast copy. */ 2074 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2075 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2076 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2077 2078 if (a->i[A->rmap->n] != b->i[B->rmap->n]) { 2079 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 2080 } 2081 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2082 } else { 2083 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2084 } 2085 PetscFunctionReturn(0); 2086 } 2087 2088 #undef __FUNCT__ 2089 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 2090 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 2091 { 2092 PetscErrorCode ierr; 2093 2094 PetscFunctionBegin; 2095 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2096 PetscFunctionReturn(0); 2097 } 2098 2099 #undef __FUNCT__ 2100 #define __FUNCT__ "MatGetArray_SeqAIJ" 2101 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2102 { 2103 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2104 PetscFunctionBegin; 2105 *array = a->a; 2106 PetscFunctionReturn(0); 2107 } 2108 2109 #undef __FUNCT__ 2110 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2111 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2112 { 2113 PetscFunctionBegin; 2114 PetscFunctionReturn(0); 2115 } 2116 2117 #undef __FUNCT__ 2118 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2119 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2120 { 2121 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2122 PetscErrorCode ierr; 2123 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2124 PetscScalar dx,*y,*xx,*w3_array; 2125 PetscScalar *vscale_array; 2126 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2127 Vec w1,w2,w3; 2128 void *fctx = coloring->fctx; 2129 PetscTruth flg = PETSC_FALSE; 2130 2131 PetscFunctionBegin; 2132 if (!coloring->w1) { 2133 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2134 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2135 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2136 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2137 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2138 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2139 } 2140 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2141 2142 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2143 ierr = PetscOptionsGetTruth(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,PETSC_NULL);CHKERRQ(ierr); 2144 if (flg) { 2145 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2146 } else { 2147 PetscTruth assembled; 2148 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2149 if (assembled) { 2150 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2151 } 2152 } 2153 2154 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2155 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2156 2157 /* 2158 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2159 coloring->F for the coarser grids from the finest 2160 */ 2161 if (coloring->F) { 2162 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2163 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2164 if (m1 != m2) { 2165 coloring->F = 0; 2166 } 2167 } 2168 2169 if (coloring->F) { 2170 w1 = coloring->F; 2171 coloring->F = 0; 2172 } else { 2173 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2174 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2175 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2176 } 2177 2178 /* 2179 Compute all the scale factors and share with other processors 2180 */ 2181 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2182 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2183 for (k=0; k<coloring->ncolors; k++) { 2184 /* 2185 Loop over each column associated with color adding the 2186 perturbation to the vector w3. 2187 */ 2188 for (l=0; l<coloring->ncolumns[k]; l++) { 2189 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2190 dx = xx[col]; 2191 if (dx == 0.0) dx = 1.0; 2192 #if !defined(PETSC_USE_COMPLEX) 2193 if (dx < umin && dx >= 0.0) dx = umin; 2194 else if (dx < 0.0 && dx > -umin) dx = -umin; 2195 #else 2196 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2197 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2198 #endif 2199 dx *= epsilon; 2200 vscale_array[col] = 1.0/dx; 2201 } 2202 } 2203 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2204 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2205 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2206 2207 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2208 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2209 2210 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2211 else vscaleforrow = coloring->columnsforrow; 2212 2213 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2214 /* 2215 Loop over each color 2216 */ 2217 for (k=0; k<coloring->ncolors; k++) { 2218 coloring->currentcolor = k; 2219 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2220 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2221 /* 2222 Loop over each column associated with color adding the 2223 perturbation to the vector w3. 2224 */ 2225 for (l=0; l<coloring->ncolumns[k]; l++) { 2226 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2227 dx = xx[col]; 2228 if (dx == 0.0) dx = 1.0; 2229 #if !defined(PETSC_USE_COMPLEX) 2230 if (dx < umin && dx >= 0.0) dx = umin; 2231 else if (dx < 0.0 && dx > -umin) dx = -umin; 2232 #else 2233 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2234 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2235 #endif 2236 dx *= epsilon; 2237 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2238 w3_array[col] += dx; 2239 } 2240 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2241 2242 /* 2243 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2244 */ 2245 2246 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2247 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2248 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2249 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2250 2251 /* 2252 Loop over rows of vector, putting results into Jacobian matrix 2253 */ 2254 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2255 for (l=0; l<coloring->nrows[k]; l++) { 2256 row = coloring->rows[k][l]; 2257 col = coloring->columnsforrow[k][l]; 2258 y[row] *= vscale_array[vscaleforrow[k][l]]; 2259 srow = row + start; 2260 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2261 } 2262 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2263 } 2264 coloring->currentcolor = k; 2265 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2266 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2267 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2268 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2269 PetscFunctionReturn(0); 2270 } 2271 2272 #undef __FUNCT__ 2273 #define __FUNCT__ "MatAXPY_SeqAIJ" 2274 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2275 { 2276 PetscErrorCode ierr; 2277 PetscInt i; 2278 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2279 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2280 2281 PetscFunctionBegin; 2282 if (str == SAME_NONZERO_PATTERN) { 2283 PetscScalar alpha = a; 2284 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2285 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2286 if (y->xtoy && y->XtoY != X) { 2287 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2288 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2289 } 2290 if (!y->xtoy) { /* get xtoy */ 2291 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2292 y->XtoY = X; 2293 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2294 } 2295 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2296 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2297 } else { 2298 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2299 } 2300 PetscFunctionReturn(0); 2301 } 2302 2303 #undef __FUNCT__ 2304 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2305 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2306 { 2307 PetscErrorCode ierr; 2308 2309 PetscFunctionBegin; 2310 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 2311 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 2312 PetscFunctionReturn(0); 2313 } 2314 2315 #undef __FUNCT__ 2316 #define __FUNCT__ "MatConjugate_SeqAIJ" 2317 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2318 { 2319 #if defined(PETSC_USE_COMPLEX) 2320 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2321 PetscInt i,nz; 2322 PetscScalar *a; 2323 2324 PetscFunctionBegin; 2325 nz = aij->nz; 2326 a = aij->a; 2327 for (i=0; i<nz; i++) { 2328 a[i] = PetscConj(a[i]); 2329 } 2330 #else 2331 PetscFunctionBegin; 2332 #endif 2333 PetscFunctionReturn(0); 2334 } 2335 2336 #undef __FUNCT__ 2337 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2338 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2339 { 2340 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2341 PetscErrorCode ierr; 2342 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2343 PetscReal atmp; 2344 PetscScalar *x; 2345 MatScalar *aa; 2346 2347 PetscFunctionBegin; 2348 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2349 aa = a->a; 2350 ai = a->i; 2351 aj = a->j; 2352 2353 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2354 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2355 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2356 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2357 for (i=0; i<m; i++) { 2358 ncols = ai[1] - ai[0]; ai++; 2359 x[i] = 0.0; 2360 for (j=0; j<ncols; j++){ 2361 atmp = PetscAbsScalar(*aa); 2362 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2363 aa++; aj++; 2364 } 2365 } 2366 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2367 PetscFunctionReturn(0); 2368 } 2369 2370 #undef __FUNCT__ 2371 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2372 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2373 { 2374 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2375 PetscErrorCode ierr; 2376 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2377 PetscScalar *x; 2378 MatScalar *aa; 2379 2380 PetscFunctionBegin; 2381 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2382 aa = a->a; 2383 ai = a->i; 2384 aj = a->j; 2385 2386 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2387 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2388 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2389 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2390 for (i=0; i<m; i++) { 2391 ncols = ai[1] - ai[0]; ai++; 2392 if (ncols == A->cmap->n) { /* row is dense */ 2393 x[i] = *aa; if (idx) idx[i] = 0; 2394 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2395 x[i] = 0.0; 2396 if (idx) { 2397 idx[i] = 0; /* in case ncols is zero */ 2398 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2399 if (aj[j] > j) { 2400 idx[i] = j; 2401 break; 2402 } 2403 } 2404 } 2405 } 2406 for (j=0; j<ncols; j++){ 2407 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2408 aa++; aj++; 2409 } 2410 } 2411 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2412 PetscFunctionReturn(0); 2413 } 2414 2415 #undef __FUNCT__ 2416 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2417 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2418 { 2419 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2420 PetscErrorCode ierr; 2421 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2422 PetscReal atmp; 2423 PetscScalar *x; 2424 MatScalar *aa; 2425 2426 PetscFunctionBegin; 2427 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2428 aa = a->a; 2429 ai = a->i; 2430 aj = a->j; 2431 2432 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2433 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2434 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2435 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2436 for (i=0; i<m; i++) { 2437 ncols = ai[1] - ai[0]; ai++; 2438 if (ncols) { 2439 /* Get first nonzero */ 2440 for(j = 0; j < ncols; j++) { 2441 atmp = PetscAbsScalar(aa[j]); 2442 if (atmp > 1.0e-12) {x[i] = atmp; if (idx) idx[i] = aj[j]; break;} 2443 } 2444 if (j == ncols) {x[i] = *aa; if (idx) idx[i] = *aj;} 2445 } else { 2446 x[i] = 0.0; if (idx) idx[i] = 0; 2447 } 2448 for(j = 0; j < ncols; j++) { 2449 atmp = PetscAbsScalar(*aa); 2450 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2451 aa++; aj++; 2452 } 2453 } 2454 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2455 PetscFunctionReturn(0); 2456 } 2457 2458 #undef __FUNCT__ 2459 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2460 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2461 { 2462 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2463 PetscErrorCode ierr; 2464 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2465 PetscScalar *x; 2466 MatScalar *aa; 2467 2468 PetscFunctionBegin; 2469 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2470 aa = a->a; 2471 ai = a->i; 2472 aj = a->j; 2473 2474 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2475 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2476 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2477 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2478 for (i=0; i<m; i++) { 2479 ncols = ai[1] - ai[0]; ai++; 2480 if (ncols == A->cmap->n) { /* row is dense */ 2481 x[i] = *aa; if (idx) idx[i] = 0; 2482 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2483 x[i] = 0.0; 2484 if (idx) { /* find first implicit 0.0 in the row */ 2485 idx[i] = 0; /* in case ncols is zero */ 2486 for (j=0;j<ncols;j++) { 2487 if (aj[j] > j) { 2488 idx[i] = j; 2489 break; 2490 } 2491 } 2492 } 2493 } 2494 for (j=0; j<ncols; j++){ 2495 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2496 aa++; aj++; 2497 } 2498 } 2499 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2500 PetscFunctionReturn(0); 2501 } 2502 extern PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 2503 /* -------------------------------------------------------------------*/ 2504 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2505 MatGetRow_SeqAIJ, 2506 MatRestoreRow_SeqAIJ, 2507 MatMult_SeqAIJ, 2508 /* 4*/ MatMultAdd_SeqAIJ, 2509 MatMultTranspose_SeqAIJ, 2510 MatMultTransposeAdd_SeqAIJ, 2511 0, 2512 0, 2513 0, 2514 /*10*/ 0, 2515 MatLUFactor_SeqAIJ, 2516 0, 2517 MatSOR_SeqAIJ, 2518 MatTranspose_SeqAIJ, 2519 /*15*/ MatGetInfo_SeqAIJ, 2520 MatEqual_SeqAIJ, 2521 MatGetDiagonal_SeqAIJ, 2522 MatDiagonalScale_SeqAIJ, 2523 MatNorm_SeqAIJ, 2524 /*20*/ 0, 2525 MatAssemblyEnd_SeqAIJ, 2526 MatSetOption_SeqAIJ, 2527 MatZeroEntries_SeqAIJ, 2528 /*24*/ MatZeroRows_SeqAIJ, 2529 0, 2530 0, 2531 0, 2532 0, 2533 /*29*/ MatSetUpPreallocation_SeqAIJ, 2534 0, 2535 0, 2536 MatGetArray_SeqAIJ, 2537 MatRestoreArray_SeqAIJ, 2538 /*34*/ MatDuplicate_SeqAIJ, 2539 0, 2540 0, 2541 MatILUFactor_SeqAIJ, 2542 0, 2543 /*39*/ MatAXPY_SeqAIJ, 2544 MatGetSubMatrices_SeqAIJ, 2545 MatIncreaseOverlap_SeqAIJ, 2546 MatGetValues_SeqAIJ, 2547 MatCopy_SeqAIJ, 2548 /*44*/ MatGetRowMax_SeqAIJ, 2549 MatScale_SeqAIJ, 2550 0, 2551 MatDiagonalSet_SeqAIJ, 2552 0, 2553 /*49*/ MatSetBlockSize_SeqAIJ, 2554 MatGetRowIJ_SeqAIJ, 2555 MatRestoreRowIJ_SeqAIJ, 2556 MatGetColumnIJ_SeqAIJ, 2557 MatRestoreColumnIJ_SeqAIJ, 2558 /*54*/ MatFDColoringCreate_SeqAIJ, 2559 0, 2560 0, 2561 MatPermute_SeqAIJ, 2562 0, 2563 /*59*/ 0, 2564 MatDestroy_SeqAIJ, 2565 MatView_SeqAIJ, 2566 0, 2567 0, 2568 /*64*/ 0, 2569 0, 2570 0, 2571 0, 2572 0, 2573 /*69*/ MatGetRowMaxAbs_SeqAIJ, 2574 MatGetRowMinAbs_SeqAIJ, 2575 0, 2576 MatSetColoring_SeqAIJ, 2577 #if defined(PETSC_HAVE_ADIC) 2578 MatSetValuesAdic_SeqAIJ, 2579 #else 2580 0, 2581 #endif 2582 /*74*/ MatSetValuesAdifor_SeqAIJ, 2583 MatFDColoringApply_AIJ, 2584 0, 2585 0, 2586 0, 2587 /*79*/ 0, 2588 0, 2589 0, 2590 0, 2591 MatLoad_SeqAIJ, 2592 /*84*/ MatIsSymmetric_SeqAIJ, 2593 MatIsHermitian_SeqAIJ, 2594 0, 2595 0, 2596 0, 2597 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 2598 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2599 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2600 MatPtAP_Basic, 2601 MatPtAPSymbolic_SeqAIJ, 2602 /*94*/ MatPtAPNumeric_SeqAIJ, 2603 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2604 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2605 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2606 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2607 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 2608 0, 2609 0, 2610 MatConjugate_SeqAIJ, 2611 0, 2612 /*104*/MatSetValuesRow_SeqAIJ, 2613 MatRealPart_SeqAIJ, 2614 MatImaginaryPart_SeqAIJ, 2615 0, 2616 0, 2617 /*109*/0, 2618 0, 2619 MatGetRowMin_SeqAIJ, 2620 0, 2621 MatMissingDiagonal_SeqAIJ, 2622 /*114*/0, 2623 0, 2624 0, 2625 0, 2626 0, 2627 /*119*/0, 2628 0, 2629 0, 2630 0, 2631 MatLoadnew_SeqAIJ 2632 }; 2633 2634 EXTERN_C_BEGIN 2635 #undef __FUNCT__ 2636 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2637 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2638 { 2639 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2640 PetscInt i,nz,n; 2641 2642 PetscFunctionBegin; 2643 2644 nz = aij->maxnz; 2645 n = mat->rmap->n; 2646 for (i=0; i<nz; i++) { 2647 aij->j[i] = indices[i]; 2648 } 2649 aij->nz = nz; 2650 for (i=0; i<n; i++) { 2651 aij->ilen[i] = aij->imax[i]; 2652 } 2653 2654 PetscFunctionReturn(0); 2655 } 2656 EXTERN_C_END 2657 2658 #undef __FUNCT__ 2659 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2660 /*@ 2661 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2662 in the matrix. 2663 2664 Input Parameters: 2665 + mat - the SeqAIJ matrix 2666 - indices - the column indices 2667 2668 Level: advanced 2669 2670 Notes: 2671 This can be called if you have precomputed the nonzero structure of the 2672 matrix and want to provide it to the matrix object to improve the performance 2673 of the MatSetValues() operation. 2674 2675 You MUST have set the correct numbers of nonzeros per row in the call to 2676 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2677 2678 MUST be called before any calls to MatSetValues(); 2679 2680 The indices should start with zero, not one. 2681 2682 @*/ 2683 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2684 { 2685 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2686 2687 PetscFunctionBegin; 2688 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2689 PetscValidPointer(indices,2); 2690 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2691 if (f) { 2692 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2693 } else { 2694 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2695 } 2696 PetscFunctionReturn(0); 2697 } 2698 2699 /* ----------------------------------------------------------------------------------------*/ 2700 2701 EXTERN_C_BEGIN 2702 #undef __FUNCT__ 2703 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2704 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2705 { 2706 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2707 PetscErrorCode ierr; 2708 size_t nz = aij->i[mat->rmap->n]; 2709 2710 PetscFunctionBegin; 2711 if (aij->nonew != 1) { 2712 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2713 } 2714 2715 /* allocate space for values if not already there */ 2716 if (!aij->saved_values) { 2717 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2718 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2719 } 2720 2721 /* copy values over */ 2722 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2723 PetscFunctionReturn(0); 2724 } 2725 EXTERN_C_END 2726 2727 #undef __FUNCT__ 2728 #define __FUNCT__ "MatStoreValues" 2729 /*@ 2730 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2731 example, reuse of the linear part of a Jacobian, while recomputing the 2732 nonlinear portion. 2733 2734 Collect on Mat 2735 2736 Input Parameters: 2737 . mat - the matrix (currently only AIJ matrices support this option) 2738 2739 Level: advanced 2740 2741 Common Usage, with SNESSolve(): 2742 $ Create Jacobian matrix 2743 $ Set linear terms into matrix 2744 $ Apply boundary conditions to matrix, at this time matrix must have 2745 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2746 $ boundary conditions again will not change the nonzero structure 2747 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2748 $ ierr = MatStoreValues(mat); 2749 $ Call SNESSetJacobian() with matrix 2750 $ In your Jacobian routine 2751 $ ierr = MatRetrieveValues(mat); 2752 $ Set nonlinear terms in matrix 2753 2754 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2755 $ // build linear portion of Jacobian 2756 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2757 $ ierr = MatStoreValues(mat); 2758 $ loop over nonlinear iterations 2759 $ ierr = MatRetrieveValues(mat); 2760 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2761 $ // call MatAssemblyBegin/End() on matrix 2762 $ Solve linear system with Jacobian 2763 $ endloop 2764 2765 Notes: 2766 Matrix must already be assemblied before calling this routine 2767 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2768 calling this routine. 2769 2770 When this is called multiple times it overwrites the previous set of stored values 2771 and does not allocated additional space. 2772 2773 .seealso: MatRetrieveValues() 2774 2775 @*/ 2776 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2777 { 2778 PetscErrorCode ierr,(*f)(Mat); 2779 2780 PetscFunctionBegin; 2781 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2782 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2783 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2784 2785 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2786 if (f) { 2787 ierr = (*f)(mat);CHKERRQ(ierr); 2788 } else { 2789 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2790 } 2791 PetscFunctionReturn(0); 2792 } 2793 2794 EXTERN_C_BEGIN 2795 #undef __FUNCT__ 2796 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2797 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2798 { 2799 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2800 PetscErrorCode ierr; 2801 PetscInt nz = aij->i[mat->rmap->n]; 2802 2803 PetscFunctionBegin; 2804 if (aij->nonew != 1) { 2805 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2806 } 2807 if (!aij->saved_values) { 2808 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2809 } 2810 /* copy values over */ 2811 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2812 PetscFunctionReturn(0); 2813 } 2814 EXTERN_C_END 2815 2816 #undef __FUNCT__ 2817 #define __FUNCT__ "MatRetrieveValues" 2818 /*@ 2819 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2820 example, reuse of the linear part of a Jacobian, while recomputing the 2821 nonlinear portion. 2822 2823 Collect on Mat 2824 2825 Input Parameters: 2826 . mat - the matrix (currently on AIJ matrices support this option) 2827 2828 Level: advanced 2829 2830 .seealso: MatStoreValues() 2831 2832 @*/ 2833 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2834 { 2835 PetscErrorCode ierr,(*f)(Mat); 2836 2837 PetscFunctionBegin; 2838 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2839 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2840 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2841 2842 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2843 if (f) { 2844 ierr = (*f)(mat);CHKERRQ(ierr); 2845 } else { 2846 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2847 } 2848 PetscFunctionReturn(0); 2849 } 2850 2851 2852 /* --------------------------------------------------------------------------------*/ 2853 #undef __FUNCT__ 2854 #define __FUNCT__ "MatCreateSeqAIJ" 2855 /*@C 2856 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2857 (the default parallel PETSc format). For good matrix assembly performance 2858 the user should preallocate the matrix storage by setting the parameter nz 2859 (or the array nnz). By setting these parameters accurately, performance 2860 during matrix assembly can be increased by more than a factor of 50. 2861 2862 Collective on MPI_Comm 2863 2864 Input Parameters: 2865 + comm - MPI communicator, set to PETSC_COMM_SELF 2866 . m - number of rows 2867 . n - number of columns 2868 . nz - number of nonzeros per row (same for all rows) 2869 - nnz - array containing the number of nonzeros in the various rows 2870 (possibly different for each row) or PETSC_NULL 2871 2872 Output Parameter: 2873 . A - the matrix 2874 2875 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2876 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2877 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2878 2879 Notes: 2880 If nnz is given then nz is ignored 2881 2882 The AIJ format (also called the Yale sparse matrix format or 2883 compressed row storage), is fully compatible with standard Fortran 77 2884 storage. That is, the stored row and column indices can begin at 2885 either one (as in Fortran) or zero. See the users' manual for details. 2886 2887 Specify the preallocated storage with either nz or nnz (not both). 2888 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2889 allocation. For large problems you MUST preallocate memory or you 2890 will get TERRIBLE performance, see the users' manual chapter on matrices. 2891 2892 By default, this format uses inodes (identical nodes) when possible, to 2893 improve numerical efficiency of matrix-vector products and solves. We 2894 search for consecutive rows with the same nonzero structure, thereby 2895 reusing matrix information to achieve increased efficiency. 2896 2897 Options Database Keys: 2898 + -mat_no_inode - Do not use inodes 2899 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2900 - -mat_aij_oneindex - Internally use indexing starting at 1 2901 rather than 0. Note that when calling MatSetValues(), 2902 the user still MUST index entries starting at 0! 2903 2904 Level: intermediate 2905 2906 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2907 2908 @*/ 2909 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2910 { 2911 PetscErrorCode ierr; 2912 2913 PetscFunctionBegin; 2914 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2915 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2916 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2917 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2918 PetscFunctionReturn(0); 2919 } 2920 2921 #undef __FUNCT__ 2922 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2923 /*@C 2924 MatSeqAIJSetPreallocation - For good matrix assembly performance 2925 the user should preallocate the matrix storage by setting the parameter nz 2926 (or the array nnz). By setting these parameters accurately, performance 2927 during matrix assembly can be increased by more than a factor of 50. 2928 2929 Collective on MPI_Comm 2930 2931 Input Parameters: 2932 + B - The matrix-free 2933 . nz - number of nonzeros per row (same for all rows) 2934 - nnz - array containing the number of nonzeros in the various rows 2935 (possibly different for each row) or PETSC_NULL 2936 2937 Notes: 2938 If nnz is given then nz is ignored 2939 2940 The AIJ format (also called the Yale sparse matrix format or 2941 compressed row storage), is fully compatible with standard Fortran 77 2942 storage. That is, the stored row and column indices can begin at 2943 either one (as in Fortran) or zero. See the users' manual for details. 2944 2945 Specify the preallocated storage with either nz or nnz (not both). 2946 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2947 allocation. For large problems you MUST preallocate memory or you 2948 will get TERRIBLE performance, see the users' manual chapter on matrices. 2949 2950 You can call MatGetInfo() to get information on how effective the preallocation was; 2951 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2952 You can also run with the option -info and look for messages with the string 2953 malloc in them to see if additional memory allocation was needed. 2954 2955 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2956 entries or columns indices 2957 2958 By default, this format uses inodes (identical nodes) when possible, to 2959 improve numerical efficiency of matrix-vector products and solves. We 2960 search for consecutive rows with the same nonzero structure, thereby 2961 reusing matrix information to achieve increased efficiency. 2962 2963 Options Database Keys: 2964 + -mat_no_inode - Do not use inodes 2965 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2966 - -mat_aij_oneindex - Internally use indexing starting at 1 2967 rather than 0. Note that when calling MatSetValues(), 2968 the user still MUST index entries starting at 0! 2969 2970 Level: intermediate 2971 2972 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 2973 2974 @*/ 2975 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2976 { 2977 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2978 2979 PetscFunctionBegin; 2980 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2981 if (f) { 2982 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2983 } 2984 PetscFunctionReturn(0); 2985 } 2986 2987 EXTERN_C_BEGIN 2988 #undef __FUNCT__ 2989 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2990 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2991 { 2992 Mat_SeqAIJ *b; 2993 PetscTruth skipallocation = PETSC_FALSE; 2994 PetscErrorCode ierr; 2995 PetscInt i; 2996 2997 PetscFunctionBegin; 2998 2999 if (nz == MAT_SKIP_ALLOCATION) { 3000 skipallocation = PETSC_TRUE; 3001 nz = 0; 3002 } 3003 3004 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 3005 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 3006 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3007 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3008 3009 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3010 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3011 if (nnz) { 3012 for (i=0; i<B->rmap->n; i++) { 3013 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]); 3014 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); 3015 } 3016 } 3017 3018 B->preallocated = PETSC_TRUE; 3019 b = (Mat_SeqAIJ*)B->data; 3020 3021 if (!skipallocation) { 3022 if (!b->imax) { 3023 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3024 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3025 } 3026 if (!nnz) { 3027 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3028 else if (nz <= 0) nz = 1; 3029 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3030 nz = nz*B->rmap->n; 3031 } else { 3032 nz = 0; 3033 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3034 } 3035 /* b->ilen will count nonzeros in each row so far. */ 3036 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3037 3038 /* allocate the matrix space */ 3039 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3040 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3041 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3042 b->i[0] = 0; 3043 for (i=1; i<B->rmap->n+1; i++) { 3044 b->i[i] = b->i[i-1] + b->imax[i-1]; 3045 } 3046 b->singlemalloc = PETSC_TRUE; 3047 b->free_a = PETSC_TRUE; 3048 b->free_ij = PETSC_TRUE; 3049 } else { 3050 b->free_a = PETSC_FALSE; 3051 b->free_ij = PETSC_FALSE; 3052 } 3053 3054 b->nz = 0; 3055 b->maxnz = nz; 3056 B->info.nz_unneeded = (double)b->maxnz; 3057 PetscFunctionReturn(0); 3058 } 3059 EXTERN_C_END 3060 3061 #undef __FUNCT__ 3062 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3063 /*@ 3064 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3065 3066 Input Parameters: 3067 + B - the matrix 3068 . i - the indices into j for the start of each row (starts with zero) 3069 . j - the column indices for each row (starts with zero) these must be sorted for each row 3070 - v - optional values in the matrix 3071 3072 Level: developer 3073 3074 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3075 3076 .keywords: matrix, aij, compressed row, sparse, sequential 3077 3078 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3079 @*/ 3080 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3081 { 3082 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 3083 PetscErrorCode ierr; 3084 3085 PetscFunctionBegin; 3086 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3087 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 3088 if (f) { 3089 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 3090 } 3091 PetscFunctionReturn(0); 3092 } 3093 3094 EXTERN_C_BEGIN 3095 #undef __FUNCT__ 3096 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3097 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3098 { 3099 PetscInt i; 3100 PetscInt m,n; 3101 PetscInt nz; 3102 PetscInt *nnz, nz_max = 0; 3103 PetscScalar *values; 3104 PetscErrorCode ierr; 3105 3106 PetscFunctionBegin; 3107 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3108 3109 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3110 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3111 for(i = 0; i < m; i++) { 3112 nz = Ii[i+1]- Ii[i]; 3113 nz_max = PetscMax(nz_max, nz); 3114 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3115 nnz[i] = nz; 3116 } 3117 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3118 ierr = PetscFree(nnz);CHKERRQ(ierr); 3119 3120 if (v) { 3121 values = (PetscScalar*) v; 3122 } else { 3123 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3124 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3125 } 3126 3127 for(i = 0; i < m; i++) { 3128 nz = Ii[i+1] - Ii[i]; 3129 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3130 } 3131 3132 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3133 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3134 3135 if (!v) { 3136 ierr = PetscFree(values);CHKERRQ(ierr); 3137 } 3138 PetscFunctionReturn(0); 3139 } 3140 EXTERN_C_END 3141 3142 #include "../src/mat/impls/dense/seq/dense.h" 3143 #include "private/petscaxpy.h" 3144 3145 #undef __FUNCT__ 3146 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3147 /* 3148 Computes (B'*A')' since computing B*A directly is untenable 3149 3150 n p p 3151 ( ) ( ) ( ) 3152 m ( A ) * n ( B ) = m ( C ) 3153 ( ) ( ) ( ) 3154 3155 */ 3156 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3157 { 3158 PetscErrorCode ierr; 3159 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3160 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3161 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3162 PetscInt i,n,m,q,p; 3163 const PetscInt *ii,*idx; 3164 const PetscScalar *b,*a,*a_q; 3165 PetscScalar *c,*c_q; 3166 3167 PetscFunctionBegin; 3168 m = A->rmap->n; 3169 n = A->cmap->n; 3170 p = B->cmap->n; 3171 a = sub_a->v; 3172 b = sub_b->a; 3173 c = sub_c->v; 3174 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3175 3176 ii = sub_b->i; 3177 idx = sub_b->j; 3178 for (i=0; i<n; i++) { 3179 q = ii[i+1] - ii[i]; 3180 while (q-->0) { 3181 c_q = c + m*(*idx); 3182 a_q = a + m*i; 3183 PetscAXPY(c_q,*b,a_q,m); 3184 idx++; 3185 b++; 3186 } 3187 } 3188 PetscFunctionReturn(0); 3189 } 3190 3191 #undef __FUNCT__ 3192 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3193 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3194 { 3195 PetscErrorCode ierr; 3196 PetscInt m=A->rmap->n,n=B->cmap->n; 3197 Mat Cmat; 3198 3199 PetscFunctionBegin; 3200 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); 3201 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3202 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3203 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3204 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3205 Cmat->assembled = PETSC_TRUE; 3206 *C = Cmat; 3207 PetscFunctionReturn(0); 3208 } 3209 3210 /* ----------------------------------------------------------------*/ 3211 #undef __FUNCT__ 3212 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3213 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3214 { 3215 PetscErrorCode ierr; 3216 3217 PetscFunctionBegin; 3218 if (scall == MAT_INITIAL_MATRIX){ 3219 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3220 } 3221 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3222 PetscFunctionReturn(0); 3223 } 3224 3225 3226 /*MC 3227 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3228 based on compressed sparse row format. 3229 3230 Options Database Keys: 3231 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3232 3233 Level: beginner 3234 3235 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3236 M*/ 3237 3238 EXTERN_C_BEGIN 3239 #if defined(PETSC_HAVE_PASTIX) 3240 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3241 #endif 3242 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3243 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3244 #endif 3245 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 3246 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3247 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3248 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscTruth *); 3249 #if defined(PETSC_HAVE_MUMPS) 3250 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3251 #endif 3252 #if defined(PETSC_HAVE_SUPERLU) 3253 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3254 #endif 3255 #if defined(PETSC_HAVE_SUPERLU_DIST) 3256 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3257 #endif 3258 #if defined(PETSC_HAVE_SPOOLES) 3259 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3260 #endif 3261 #if defined(PETSC_HAVE_UMFPACK) 3262 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3263 #endif 3264 #if defined(PETSC_HAVE_CHOLMOD) 3265 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3266 #endif 3267 #if defined(PETSC_HAVE_LUSOL) 3268 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3269 #endif 3270 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3271 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3272 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEnginePut_SeqAIJ(PetscObject,void*); 3273 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEngineGet_SeqAIJ(PetscObject,void*); 3274 #endif 3275 EXTERN_C_END 3276 3277 3278 EXTERN_C_BEGIN 3279 #undef __FUNCT__ 3280 #define __FUNCT__ "MatCreate_SeqAIJ" 3281 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3282 { 3283 Mat_SeqAIJ *b; 3284 PetscErrorCode ierr; 3285 PetscMPIInt size; 3286 3287 PetscFunctionBegin; 3288 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3289 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3290 3291 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3292 B->data = (void*)b; 3293 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3294 B->mapping = 0; 3295 b->row = 0; 3296 b->col = 0; 3297 b->icol = 0; 3298 b->reallocs = 0; 3299 b->ignorezeroentries = PETSC_FALSE; 3300 b->roworiented = PETSC_TRUE; 3301 b->nonew = 0; 3302 b->diag = 0; 3303 b->solve_work = 0; 3304 B->spptr = 0; 3305 b->saved_values = 0; 3306 b->idiag = 0; 3307 b->mdiag = 0; 3308 b->ssor_work = 0; 3309 b->omega = 1.0; 3310 b->fshift = 0.0; 3311 b->idiagvalid = PETSC_FALSE; 3312 b->keepnonzeropattern = PETSC_FALSE; 3313 b->xtoy = 0; 3314 b->XtoY = 0; 3315 b->compressedrow.use = PETSC_FALSE; 3316 b->compressedrow.nrows = B->rmap->n; 3317 b->compressedrow.i = PETSC_NULL; 3318 b->compressedrow.rindex = PETSC_NULL; 3319 b->compressedrow.checked = PETSC_FALSE; 3320 B->same_nonzero = PETSC_FALSE; 3321 3322 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3323 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3324 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C", 3325 "MatGetFactor_seqaij_matlab", 3326 MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3327 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3328 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3329 #endif 3330 #if defined(PETSC_HAVE_PASTIX) 3331 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 3332 "MatGetFactor_seqaij_pastix", 3333 MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3334 #endif 3335 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3336 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C", 3337 "MatGetFactor_seqaij_essl", 3338 MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3339 #endif 3340 #if defined(PETSC_HAVE_SUPERLU) 3341 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C", 3342 "MatGetFactor_seqaij_superlu", 3343 MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3344 #endif 3345 #if defined(PETSC_HAVE_SUPERLU_DIST) 3346 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C", 3347 "MatGetFactor_seqaij_superlu_dist", 3348 MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3349 #endif 3350 #if defined(PETSC_HAVE_SPOOLES) 3351 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 3352 "MatGetFactor_seqaij_spooles", 3353 MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3354 #endif 3355 #if defined(PETSC_HAVE_MUMPS) 3356 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 3357 "MatGetFactor_aij_mumps", 3358 MatGetFactor_aij_mumps);CHKERRQ(ierr); 3359 #endif 3360 #if defined(PETSC_HAVE_UMFPACK) 3361 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C", 3362 "MatGetFactor_seqaij_umfpack", 3363 MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3364 #endif 3365 #if defined(PETSC_HAVE_CHOLMOD) 3366 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 3367 "MatGetFactor_seqaij_cholmod", 3368 MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3369 #endif 3370 #if defined(PETSC_HAVE_LUSOL) 3371 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C", 3372 "MatGetFactor_seqaij_lusol", 3373 MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3374 #endif 3375 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 3376 "MatGetFactor_seqaij_petsc", 3377 MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3378 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 3379 "MatGetFactorAvailable_seqaij_petsc", 3380 MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3381 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C", 3382 "MatGetFactor_seqaij_bas", 3383 MatGetFactor_seqaij_bas); 3384 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3385 "MatSeqAIJSetColumnIndices_SeqAIJ", 3386 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3387 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3388 "MatStoreValues_SeqAIJ", 3389 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3390 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3391 "MatRetrieveValues_SeqAIJ", 3392 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3393 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3394 "MatConvert_SeqAIJ_SeqSBAIJ", 3395 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3396 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3397 "MatConvert_SeqAIJ_SeqBAIJ", 3398 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3399 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3400 "MatConvert_SeqAIJ_SeqCSRPERM", 3401 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3402 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3403 "MatConvert_SeqAIJ_SeqCRL", 3404 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3405 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3406 "MatIsTranspose_SeqAIJ", 3407 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3408 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3409 "MatIsHermitianTranspose_SeqAIJ", 3410 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3411 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3412 "MatSeqAIJSetPreallocation_SeqAIJ", 3413 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3414 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3415 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3416 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3417 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3418 "MatReorderForNonzeroDiagonal_SeqAIJ", 3419 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3420 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C", 3421 "MatMatMult_SeqDense_SeqAIJ", 3422 MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3423 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C", 3424 "MatMatMultSymbolic_SeqDense_SeqAIJ", 3425 MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3426 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C", 3427 "MatMatMultNumeric_SeqDense_SeqAIJ", 3428 MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3429 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3430 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3431 PetscFunctionReturn(0); 3432 } 3433 EXTERN_C_END 3434 3435 #undef __FUNCT__ 3436 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3437 /* 3438 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3439 */ 3440 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscTruth mallocmatspace) 3441 { 3442 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3443 PetscErrorCode ierr; 3444 PetscInt i,m = A->rmap->n; 3445 3446 PetscFunctionBegin; 3447 c = (Mat_SeqAIJ*)C->data; 3448 3449 C->factortype = A->factortype; 3450 c->row = 0; 3451 c->col = 0; 3452 c->icol = 0; 3453 c->reallocs = 0; 3454 3455 C->assembled = PETSC_TRUE; 3456 3457 ierr = PetscLayoutSetBlockSize(C->rmap,1);CHKERRQ(ierr); 3458 ierr = PetscLayoutSetBlockSize(C->cmap,1);CHKERRQ(ierr); 3459 ierr = PetscLayoutSetUp(C->rmap);CHKERRQ(ierr); 3460 ierr = PetscLayoutSetUp(C->cmap);CHKERRQ(ierr); 3461 3462 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3463 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3464 for (i=0; i<m; i++) { 3465 c->imax[i] = a->imax[i]; 3466 c->ilen[i] = a->ilen[i]; 3467 } 3468 3469 /* allocate the matrix space */ 3470 if (mallocmatspace){ 3471 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3472 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3473 c->singlemalloc = PETSC_TRUE; 3474 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3475 if (m > 0) { 3476 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3477 if (cpvalues == MAT_COPY_VALUES) { 3478 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3479 } else { 3480 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3481 } 3482 } 3483 } 3484 3485 c->ignorezeroentries = a->ignorezeroentries; 3486 c->roworiented = a->roworiented; 3487 c->nonew = a->nonew; 3488 if (a->diag) { 3489 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3490 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3491 for (i=0; i<m; i++) { 3492 c->diag[i] = a->diag[i]; 3493 } 3494 } else c->diag = 0; 3495 c->solve_work = 0; 3496 c->saved_values = 0; 3497 c->idiag = 0; 3498 c->ssor_work = 0; 3499 c->keepnonzeropattern = a->keepnonzeropattern; 3500 c->free_a = PETSC_TRUE; 3501 c->free_ij = PETSC_TRUE; 3502 c->xtoy = 0; 3503 c->XtoY = 0; 3504 3505 c->nz = a->nz; 3506 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3507 C->preallocated = PETSC_TRUE; 3508 3509 c->compressedrow.use = a->compressedrow.use; 3510 c->compressedrow.nrows = a->compressedrow.nrows; 3511 c->compressedrow.checked = a->compressedrow.checked; 3512 if (a->compressedrow.checked && a->compressedrow.use){ 3513 i = a->compressedrow.nrows; 3514 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3515 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3516 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3517 } else { 3518 c->compressedrow.use = PETSC_FALSE; 3519 c->compressedrow.i = PETSC_NULL; 3520 c->compressedrow.rindex = PETSC_NULL; 3521 } 3522 C->same_nonzero = A->same_nonzero; 3523 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3524 3525 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3526 PetscFunctionReturn(0); 3527 } 3528 3529 #undef __FUNCT__ 3530 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3531 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3532 { 3533 PetscErrorCode ierr; 3534 3535 PetscFunctionBegin; 3536 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3537 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3538 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3539 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3540 PetscFunctionReturn(0); 3541 } 3542 3543 #undef __FUNCT__ 3544 #define __FUNCT__ "MatLoad_SeqAIJ" 3545 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, const MatType type,Mat *A) 3546 { 3547 Mat_SeqAIJ *a; 3548 Mat B; 3549 PetscErrorCode ierr; 3550 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 3551 int fd; 3552 PetscMPIInt size; 3553 MPI_Comm comm; 3554 3555 PetscFunctionBegin; 3556 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3557 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3558 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3559 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3560 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3561 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3562 M = header[1]; N = header[2]; nz = header[3]; 3563 3564 if (nz < 0) { 3565 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3566 } 3567 3568 /* read in row lengths */ 3569 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3570 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3571 3572 /* check if sum of rowlengths is same as nz */ 3573 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3574 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); 3575 3576 /* create our matrix */ 3577 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 3578 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3579 ierr = MatSetType(B,type);CHKERRQ(ierr); 3580 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 3581 a = (Mat_SeqAIJ*)B->data; 3582 3583 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3584 3585 /* read in nonzero values */ 3586 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3587 3588 /* set matrix "i" values */ 3589 a->i[0] = 0; 3590 for (i=1; i<= M; i++) { 3591 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3592 a->ilen[i-1] = rowlengths[i-1]; 3593 } 3594 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3595 3596 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3597 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3598 *A = B; 3599 PetscFunctionReturn(0); 3600 } 3601 3602 #undef __FUNCT__ 3603 #define __FUNCT__ "MatLoadnew_SeqAIJ" 3604 PetscErrorCode MatLoadnew_SeqAIJ(PetscViewer viewer,Mat newMat) 3605 { 3606 Mat_SeqAIJ *a; 3607 PetscErrorCode ierr; 3608 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 3609 int fd; 3610 PetscMPIInt size; 3611 MPI_Comm comm; 3612 3613 PetscFunctionBegin; 3614 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3615 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3616 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3617 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3618 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3619 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3620 M = header[1]; N = header[2]; nz = header[3]; 3621 3622 if (nz < 0) { 3623 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3624 } 3625 3626 /* read in row lengths */ 3627 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3628 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3629 3630 /* check if sum of rowlengths is same as nz */ 3631 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3632 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); 3633 3634 /* set global size if not set already*/ 3635 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 3636 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3637 } else { 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