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