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