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