Bug Summary

File:ksp/ksp/impls/cg/mpcg/mpcg.c
Warning:line 141, column 14
The left operand of '==' is a garbage value

Annotated Source Code

[?] Use j/k keys for keyboard navigation

1/*
2 This file defines the multipreconditioned conjugate gradient solver.
3
4 It is adapted to SPD matrices (or even more general thanks to full orthogonalization),
5 and any preconditioner which provides the PCApplyMP method
6*/
7#include <petsc/private/kspimpl.h>
8#include <../src/mat/impls/dense/seq/dense.h>
9#include <petscblaslapack.h>
10#include <petscmat.h>
11
12
13typedef struct {
14 PetscScalar eigtrunc; /*Threshold for eigenvalue-based pseudo inversion*/
15} KSP_MPCG;
16
17
18/*
19 KSPSetUp_MPCG - Sets up the workspace needed by the MPCG method.
20
21 This is called once, usually automatically by KSPSolve() or KSPSetUp()
22 but can be called directly by KSPSetUp()
23*/
24
25
26
27static PetscErrorCode KSPSetUp_MPCG(KSP ksp)
28{
29 PetscErrorCode ierr;
30
31 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
; petscstack->line[petscstack->currentsize] = 31; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);
32 /* get work vectors needed by MPCG */
33 ierr = KSPSetWorkVecs(ksp, 1);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),33,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
34 PetscFunctionReturn(0)do { do { ; if (petscstack && petscstack->currentsize
> 0) { petscstack->currentsize--; petscstack->function
[petscstack->currentsize] = 0; petscstack->file[petscstack
->currentsize] = 0; petscstack->line[petscstack->currentsize
] = 0; petscstack->petscroutine[petscstack->currentsize
] = PETSC_FALSE; } if (petscstack) { petscstack->hotdepth =
(((petscstack->hotdepth-1)<(0)) ? (0) : (petscstack->
hotdepth-1)); } ; } while (0); return(0);} while (0);
35}
36
37/*
38 KSPSolve_MPCG - This routine actually applies the multipreconditioned conjugate gradient method
39
40 Input Parameter:
41 . ksp - the Krylov space object that was set to use multipreconditioned conjugate gradient, by, for
42 example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPMPCG);
43*/
44static PetscErrorCode KSPSolve_MPCG(KSP ksp)
45{
46 PetscErrorCode ierr;
47 PetscInt i, comm_size;
48 PetscInt n_sd, n_coarse;
49 PetscInt VecLocSize, VecGloSize;
50 PetscInt nkk,ngood,ii,jj;
51 PetscBLASInt info, n;
52 PetscScalar beta;
1
'beta' declared without an initial value
→
53 PetscScalar *thediag, *myrows;
54 PetscReal dp = 0.0;
55 Vec X, B, R;
56 Vec Vdiag, Vsubdiag;
57 Vec Beta,Beta2;
58 Mat Z, W, P, Ptemp, Ztest, Vmat;
59 Mat Gamma, Delta, Delta_R, Delta_RD;
60 Mat *Pstor, *Wstor;
61 Mat Amat, Pmat;
62 Mat_SeqDense *mat;
63 MatScalar *vals;
64 PetscBool diagonalscale;
65 IS rows,cols;
66 const PetscInt *irows,*icols;
67
68 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
; petscstack->line[petscstack->currentsize] = 68; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);
69
70 ierr = PCGetDiagonalScale(ksp->pc, &diagonalscale);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),70,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
71 if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name)return PetscError(PetscObjectComm((PetscObject)ksp),71,__func__
,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,56,PETSC_ERROR_INITIAL,"Krylov method %s does not support diagonal scaling"
,((PetscObject)ksp)->type_name);
2
←
Assuming 'diagonalscale' is 0
→
3
←
Taking false branch
→
72
73 X = ksp->vec_sol;
74 B = ksp->vec_rhs;
75 R = ksp->work[0];
76
77 ierr = PetscMalloc2(ksp->max_it, &Pstor, ksp->max_it, &Wstor)PetscMallocA(2,PETSC_FALSE,77,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,(size_t)(ksp->max_it)*sizeof(**(&Pstor)),(&Pstor)
,(size_t)(ksp->max_it)*sizeof(**(&Wstor)),(&Wstor)
);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),77,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
78 ierr = MPI_Comm_size(((PetscObject)X)->comm, &comm_size);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),78,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
79
80 /* Ugly reverse communication to get info from the multipreconditioner */
81 ierr = MatCreateSeqDense(PETSC_COMM_SELF((MPI_Comm)0x44000001), 1, 3, NULL((void*)0), &Ztest);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),81,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
82 ierr = KSP_PCApplyMP(ksp, R, Ztest);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),82,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
83 ierr = MatDenseGetArray(Ztest, &vals);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),83,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
84 n_sd = (PetscInt)vals[0];
85 n_coarse = (PetscInt)vals[2];
86 ierr = MatDenseRestoreArray(Ztest, &vals);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),86,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
87 ierr = MatDestroy(&Ztest);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),87,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
88 /* End of ugly reverse communication to get info from the multipreconditioner */
89
90 ierr = VecCreateSeq(MPI_COMM_SELF((MPI_Comm)0x44000001), n_sd + n_coarse, &Vdiag);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),90,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
91 ierr = VecGetArray(Vdiag, &thediag);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),91,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
92 ierr = VecGetLocalSize(R, &VecLocSize);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),92,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
93 ierr = VecGetSize(R, &VecGloSize);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),93,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
94 n_sd += n_coarse;
95 ierr = MatCreateAIJ(((PetscObject)X)->comm, VecLocSize, PETSC_DECIDE-1, VecGloSize, n_sd, n_sd, NULL((void*)0), n_sd, NULL((void*)0), &Z);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),95,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
96 ierr = VecCreateMPI(((PetscObject)X)->comm, PETSC_DECIDE-1, n_sd, &Beta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),96,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
97 ierr = MatCreateSeqDense(PETSC_COMM_SELF((MPI_Comm)0x44000001), n_sd, n_sd, NULL((void*)0), &Delta_RD);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),97,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
98 ierr = PCGetOperators(ksp->pc, &Amat, &Pmat);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),98,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
99
100 ksp->its = 0;
101 if (!ksp->guess_zero) {
4
←
Assuming the condition is false
→
5
←
Taking false branch
→
102 ierr = KSP_MatMult(ksp, Amat, X, R);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),102,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* r <- b - Ax */
103 ierr = VecAYPX(R, -1.0, B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),103,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
104 } else {
105 ierr = VecCopy(B, R);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),105,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* r <- b (x is 0) */
106 }
107
108 switch (ksp->normtype) {
6
←
Control jumps to 'case KSP_NORM_UNPRECONDITIONED:' at line 115
→
109 case KSP_NORM_NATURAL:
110 ierr = KSP_PCApplyMP(ksp, R, Z);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),110,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* Z <- [B1 r, B2 r ... ] */
111 ierr = MatMultTranspose(Z, R, Beta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),111,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
112 ierr = VecSum(Beta, &beta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),112,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* Beta <- Z'*r */
113 dp = PetscSqrtReal(PetscAbsScalar(beta))sqrt(PetscAbsScalar(beta)); /* beta <- r'*(sum Bi)*r */
114 break;
115 case KSP_NORM_UNPRECONDITIONED:
116 ierr = VecNorm(R, NORM_2, &dp);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),116,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* dp <- r'*r = e'*A'*A*e */
117 break;
7
←
 Execution continues on line 124
→
118 case KSP_NORM_NONE:
119 dp = 0.0;
120 default:
121 SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype])return PetscError(PetscObjectComm((PetscObject)ksp),121,__func__
,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,56,PETSC_ERROR_INITIAL,"%s",KSPNormTypes[ksp->normtype]);
122 }
123
124 ierr = KSPLogResidualHistory(ksp, dp);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),124,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
125 ierr = KSPMonitor(ksp, 0, dp);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),125,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
126 ksp->rnorm = dp;
127 ierr = (*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),127,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* test for convergence */
128 if (ksp->reason) PetscFunctionReturn(0)do { do { ; if (petscstack && petscstack->currentsize
> 0) { petscstack->currentsize--; petscstack->function
[petscstack->currentsize] = 0; petscstack->file[petscstack
->currentsize] = 0; petscstack->line[petscstack->currentsize
] = 0; petscstack->petscroutine[petscstack->currentsize
] = PETSC_FALSE; } if (petscstack) { petscstack->hotdepth =
(((petscstack->hotdepth-1)<(0)) ? (0) : (petscstack->
hotdepth-1)); } ; } while (0); return(0);} while (0);
8
←
Assuming the condition is false
→
9
←
Taking false branch
→
129
130 if (ksp->normtype != KSP_NORM_NATURAL) {
10
←
Assuming the condition is false
→
11
←
Taking false branch
→
131 ierr = KSP_PCApplyMP(ksp, R, Z);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),131,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* Z <- [B1 r, B2 r ... ] */
132 ierr = MatMultTranspose(Z, R, Beta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),132,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
133 ierr = VecSum(Beta, &beta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),133,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* Beta <- Z'*r */
134 }
135
136 ierr = MatDuplicate(Z, MAT_COPY_VALUES, &P);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),136,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* P <- Z */
137
138 i = 0;
139 do {
140 ksp->its = i + 1;
141 if (beta == 0.0) {
12
←
The left operand of '==' is a garbage value
142 ksp->reason = KSP_CONVERGED_ATOL;
143 ierr = PetscInfo(ksp, "converged due to beta = 0\n")PetscInfo_Private(__func__,ksp,"converged due to beta = 0\n");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),143,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
144 break;
145 }
146
147 ierr = KSP_MatMatMult(ksp, Amat, P, &W, MAT_INITIAL_MATRIX);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),147,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
148 ierr = MatTransposeMatMult(P, W, MAT_INITIAL_MATRIX, PETSC_DEFAULT-2, &Delta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),148,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* Delta <- P' * W */
149
150 /* A-orthonormalization of P and W wrt all previous search direction */
151 /* This version makes use of dense algebra should be improved for larger blocks */
152
153 ierr = MatCreateRedundantMatrix(Delta, comm_size, MPI_COMM_NULL((MPI_Comm)0x04000000), MAT_INITIAL_MATRIX, &Delta_R);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),153,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
154 ierr = MatDestroy(&Delta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),154,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
155 ierr = MatConvert(Delta_R, MATSEQDENSE"seqdense", MAT_REUSE_MATRIX, &Delta_RD);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),155,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
156 ierr = MatDestroy(&Delta_R);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),156,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
157
158 mat = (Mat_SeqDense *)Delta_RD->data;
159
160 ierr = PetscBLASIntCast(Delta_RD->cmap->n, &n);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),160,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
161
162 if (!mat->fwork) {
163 PetscScalar dummy;
164 mat->lfwork = -1;
165 PetscStackCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &n, mat->v, &mat->lda, thediag, &dummy, &mat->lfwork, &info))do { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = "LAPACKsyev"; petscstack->file[petscstack->currentsize
] = "/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
; petscstack->line[petscstack->currentsize] = 165; petscstack
->petscroutine[petscstack->currentsize] = PETSC_FALSE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_TRUE || petscstack->hotdepth); } ; } while (0);
dsyev_("V", "U", &n, mat->v, &mat->lda, thediag
, &dummy, &mat->lfwork, &info); do { do {PetscErrorCode
_7_ierr = PetscMallocValidate(165,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
);do {if (__builtin_expect(!!(_7_ierr),0)) return PetscError(
((MPI_Comm)0x44000001),165,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,_7_ierr,PETSC_ERROR_REPEAT," ");} while (0);} while(0); do {
; if (petscstack && petscstack->currentsize > 0
) { petscstack->currentsize--; petscstack->function[petscstack
->currentsize] = 0; petscstack->file[petscstack->currentsize
] = 0; petscstack->line[petscstack->currentsize] = 0; petscstack
->petscroutine[petscstack->currentsize] = PETSC_FALSE; }
if (petscstack) { petscstack->hotdepth = (((petscstack->
hotdepth-1)<(0)) ? (0) : (petscstack->hotdepth-1)); } ;
} while (0); } while (0); } while(0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),165,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
166
167 mat->lfwork = (PetscInt)PetscRealPart(dummy)(dummy);
168 ierr = PetscMalloc1(mat->lfwork, &mat->fwork)PetscMallocA(1,PETSC_FALSE,168,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,(size_t)(mat->lfwork)*sizeof(**(&mat->fwork)),(&
mat->fwork));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),168,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
169 ierr = PetscLogObjectMemory((PetscObject)Delta_RD, mat->lfwork * sizeof(PetscBLASInt));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),169,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
170 }
171 PetscStackCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &n, mat->v, &mat->lda, thediag, mat->fwork, &mat->lfwork, &info))do { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = "LAPACKsyev"; petscstack->file[petscstack->currentsize
] = "/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
; petscstack->line[petscstack->currentsize] = 171; petscstack
->petscroutine[petscstack->currentsize] = PETSC_FALSE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_TRUE || petscstack->hotdepth); } ; } while (0);
dsyev_("V", "U", &n, mat->v, &mat->lda, thediag
, mat->fwork, &mat->lfwork, &info); do { do {PetscErrorCode
_7_ierr = PetscMallocValidate(171,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
);do {if (__builtin_expect(!!(_7_ierr),0)) return PetscError(
((MPI_Comm)0x44000001),171,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,_7_ierr,PETSC_ERROR_REPEAT," ");} while (0);} while(0); do {
; if (petscstack && petscstack->currentsize > 0
) { petscstack->currentsize--; petscstack->function[petscstack
->currentsize] = 0; petscstack->file[petscstack->currentsize
] = 0; petscstack->line[petscstack->currentsize] = 0; petscstack
->petscroutine[petscstack->currentsize] = PETSC_FALSE; }
if (petscstack) { petscstack->hotdepth = (((petscstack->
hotdepth-1)<(0)) ? (0) : (petscstack->hotdepth-1)); } ;
} while (0); } while (0); } while(0);
172
173 if (info < 0)
174 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Bad parameter %D", (PetscInt)info - 1)return PetscError(((MPI_Comm)0x44000001),174,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,81,PETSC_ERROR_INITIAL,"Bad parameter %D",(PetscInt)info - 1
);
175
176 /* Now we analyze the eigenvalues, they are in ascending order, they also should be positive*/
177
178 nkk = 0;
179 for (jj=0; jj<n; ++jj) {
180 if (thediag[jj] < thediag[n - 1] * ((KSP_MPCG *)ksp->data)->eigtrunc) {
181 nkk = jj + 1;
182 } else {
183 thediag[jj] = 1. / sqrt(thediag[jj]);
184 }
185 }
186 ngood = n - nkk;
187
188 /* Values from nkk to n should be kept */
189
190 ierr = ISCreateStride(PetscObjectComm((PetscObject)ksp), n_sd, 0, 1, &rows);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),190,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
191 ierr = ISCreateStride(PetscObjectComm((PetscObject)ksp), ngood, 0, 1, &cols);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),191,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
192 ierr = ISGetIndices(rows,&irows);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),192,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
193 ierr = ISGetIndices(cols,&icols);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),193,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
194
195 ierr = VecCreateMPI(PetscObjectComm((PetscObject)ksp), PETSC_DECIDE-1, ngood,&Vsubdiag);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),195,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
196 ierr = VecSetValues(Vsubdiag, ngood, icols, &thediag[nkk], INSERT_VALUES);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),196,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
197 ierr = VecAssemblyBegin(Vsubdiag);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),197,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
198 ierr = VecAssemblyEnd(Vsubdiag);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),198,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
199
200 ierr = VecCreateMPI(PetscObjectComm((PetscObject)ksp), PETSC_DECIDE-1, ngood, &Beta2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),200,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
201
202 PetscMalloc1(ngood*n_sd,&myrows)PetscMallocA(1,PETSC_FALSE,202,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,(size_t)(ngood*n_sd)*sizeof(**(&myrows)),(&myrows));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),202,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
203 for (ii=0;ii<n;++ii){
204 for (jj=0;jj<ngood;++jj){
205 myrows[jj+ii*ngood]=mat->v[(nkk+jj)*n+ii];
206 }
207 }
208
209 ierr = MatCreateDense(MPI_COMM_WORLD((MPI_Comm)0x44000000), PETSC_DECIDE-1, PETSC_DECIDE-1, n, ngood, NULL((void*)0), &Vmat);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),209,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
210 ierr = MatSetValues(Vmat,n_sd,irows,ngood,icols,myrows,INSERT_VALUES);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),210,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
211 PetscFree(myrows)((*PetscTrFree)((void*)(myrows),211,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
) || ((myrows) = 0,0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),211,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
212
213 ierr = MatAssemblyBegin(Vmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),213,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
214 ierr = MatAssemblyEnd(Vmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),214,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
215 ierr = ISRestoreIndices(rows,&irows);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),215,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
216 ierr = ISRestoreIndices(cols,&icols);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),216,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
217 ierr = ISDestroy(&rows);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),217,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
218 ierr = ISDestroy(&cols);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),218,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
219
220 ierr = MatConvert(Vmat, MATAIJ"aij", MAT_INPLACE_MATRIX, &Vmat);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),220,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
221 ierr = MatMatMult(P, Vmat, MAT_INITIAL_MATRIX, PETSC_DEFAULT-2, &Pstor[i]);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),221,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
222 ierr = MatMatMult(W, Vmat, MAT_INITIAL_MATRIX, PETSC_DEFAULT-2, &Wstor[i]);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),222,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
223
224 ierr = MatDiagonalScale(Pstor[i], NULL((void*)0), Vsubdiag);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),224,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
225 ierr = MatDiagonalScale(Wstor[i], NULL((void*)0), Vsubdiag);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),225,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
226
227 ierr = MatMultTranspose(Vmat, Beta, Beta2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),227,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
228 ierr = VecPointwiseMult(Beta2, Vsubdiag, Beta2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),228,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
229 ierr = VecDestroy(&Vsubdiag);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),229,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
230 ierr = MatDestroy(&Vmat);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),230,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
231
232 ierr = MatMultAdd(Pstor[i], Beta2, X, X);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),232,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
233 ierr = VecScale(Beta2, -1.);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),233,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
234 ierr = MatMultAdd(Wstor[i], Beta2, R, R);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),234,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
235
236 ierr = MatConvert(Pstor[i], MATAIJ"aij", MAT_INPLACE_MATRIX, &Pstor[i]);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),236,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
237 ierr = MatConvert(Wstor[i], MATAIJ"aij", MAT_INPLACE_MATRIX, &Wstor[i]);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),237,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
238
239 ierr = MatDestroy(&W);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),239,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
240 ierr = MatDestroy(&P);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),240,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
241 ierr = VecDestroy(&Beta2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),241,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
242
243 ierr = KSP_PCApplyMP(ksp, R, Z);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),243,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* Z <- [B1 r, B2 r ... ] */
244 ierr = MatDuplicate(Z, MAT_COPY_VALUES, &P);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),244,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* P <- Z */
245
246 for (jj = 0; jj < i + 1; jj++) { /* Full block orthogonalization */
247 ierr = MatTransposeMatMult(Wstor[jj], P, MAT_INITIAL_MATRIX, PETSC_DEFAULT-2, &Gamma);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),247,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
248 ierr = MatMatMult(Pstor[jj], Gamma, MAT_INITIAL_MATRIX, PETSC_DEFAULT-2, &Ptemp);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),248,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
249 ierr = MatDestroy(&Gamma);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),249,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
250 ierr = MatAXPY(P, -1., Ptemp, DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),250,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
251 ierr = MatDestroy(&Ptemp);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),251,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
252 }
253
254 ierr = MatMultTranspose(P, R, Beta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),254,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
255 ierr = VecSum(Beta, &beta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),255,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* Beta <- Z'*r */
256 dp = PetscSqrtReal(PetscAbsScalar(beta))sqrt(PetscAbsScalar(beta));
257
258 if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
259 ierr = VecNorm(R, NORM_2, &dp);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),259,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0); /* dp <- r'*r */
260 }
261
262 ksp->rnorm = dp;
263 ierr = KSPLogResidualHistory(ksp, dp);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),263,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
264 ierr = KSPMonitor(ksp, i + 1, dp);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),264,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
265
266 ierr = (*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),266,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
267
268 if (ksp->reason) {
269 break;
270 }
271 i++;
272 } while (i < ksp->max_it);
273
274 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
275
276 ierr = VecRestoreArray(Vdiag, &thediag);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),276,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
277 ierr = VecDestroy(&Vdiag);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),277,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
278 ierr = VecDestroy(&Beta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),278,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
279
280 ierr = MatDestroy(&Delta_RD);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),280,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
281 ierr = MatDestroy(&Delta_R);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),281,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
282 ierr = MatDestroy(&Delta);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),282,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
283 ierr = MatDestroy(&Gamma);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),283,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
284 ierr = MatDestroy(&P);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),284,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
285 ierr = MatDestroy(&W);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),285,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
286 ierr = MatDestroy(&Z);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),286,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
287
288 i = PetscMin(i + 1, ksp->max_it)(((i + 1)<(ksp->max_it)) ? (i + 1) : (ksp->max_it));
289 for (jj = 0; jj < i; jj++) {
290 ierr = MatDestroy(&Pstor[jj]);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),290,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
291 ierr = MatDestroy(&Wstor[jj]);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),291,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
292 }
293 ierr = PetscFree2(Pstor,Wstor)PetscFreeA(2,293,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,&(Pstor),&(Wstor));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),293,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
294
295 PetscFunctionReturn(0)do { do { ; if (petscstack && petscstack->currentsize
> 0) { petscstack->currentsize--; petscstack->function
[petscstack->currentsize] = 0; petscstack->file[petscstack
->currentsize] = 0; petscstack->line[petscstack->currentsize
] = 0; petscstack->petscroutine[petscstack->currentsize
] = PETSC_FALSE; } if (petscstack) { petscstack->hotdepth =
(((petscstack->hotdepth-1)<(0)) ? (0) : (petscstack->
hotdepth-1)); } ; } while (0); return(0);} while (0);
296}
297
298/*
299 KSPSetFromOptions_MPCG - Checks the options database for options related to the
300 multipreconditioned conjugate gradient method.
301*/
302PetscErrorCode KSPSetFromOptions_MPCG(PetscOptionItems *PetscOptionsObject, KSP ksp)
303{
304 PetscErrorCode ierr;
305 PetscBool flg;
306 KSP_MPCG *cg = (KSP_MPCG *)ksp->data;
307
308 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
; petscstack->line[petscstack->currentsize] = 308; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);
309 ierr = PetscOptionsHead(PetscOptionsObject, "KSP MPCG options");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),309,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
310 cg->eigtrunc = 1.e-12;
311 ierr = PetscOptionsScalar("-ksp_mp_eigtrunc", "Truncation for pseudoinversion based on eigenvalue", "KSPMPC", cg->eigtrunc, &cg->eigtrunc, &flg)PetscOptionsScalar_Private(PetscOptionsObject,"-ksp_mp_eigtrunc"
,"Truncation for pseudoinversion based on eigenvalue","KSPMPC"
,cg->eigtrunc,&cg->eigtrunc,&flg);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),311,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
312 ierr = PetscOptionsTail()0; do {if (PetscOptionsObject->count != 1) do { do { ; if (
petscstack && petscstack->currentsize > 0) { petscstack
->currentsize--; petscstack->function[petscstack->currentsize
] = 0; petscstack->file[petscstack->currentsize] = 0; petscstack
->line[petscstack->currentsize] = 0; petscstack->petscroutine
[petscstack->currentsize] = PETSC_FALSE; } if (petscstack)
{ petscstack->hotdepth = (((petscstack->hotdepth-1)<
(0)) ? (0) : (petscstack->hotdepth-1)); } ; } while (0); return
(0);} while (0);} while(0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),312,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
313 PetscFunctionReturn(0)do { do { ; if (petscstack && petscstack->currentsize
> 0) { petscstack->currentsize--; petscstack->function
[petscstack->currentsize] = 0; petscstack->file[petscstack
->currentsize] = 0; petscstack->line[petscstack->currentsize
] = 0; petscstack->petscroutine[petscstack->currentsize
] = PETSC_FALSE; } if (petscstack) { petscstack->hotdepth =
(((petscstack->hotdepth-1)<(0)) ? (0) : (petscstack->
hotdepth-1)); } ; } while (0); return(0);} while (0);
314}
315
316/*MC
317 KSPMPCG - Multipreconditioned conjugate gradient method.
318
319 At each iteration this method searches the approximation in a large subspace which is defined by the preconditioner.
320 The preconditioner must return a Matrix instead of a Vector.
321
322 Level: intermediate
323
324 Notes: This is a simple implementation of multipreconditioned conjugate gradient.
325 It is adapted to any preconditioner providing the PCApplyMP method.
326 The computation of the pseudo inverse is based on the computation of eigenvalues, as proposed in [Molina and Roux. doi:10.1002/nme.6024].
327 This implementation is based on dense algebra for the pseudo-inversion, it is general but can not scale well with the number of columns.
328 The method proved to be efficient compared to cg for poorly conditioned system solved with FETI method [Bovet et al. doi:10.1016/j.crma.2017.01.010].
329 Depending on the PC context, optimization could be considered:
330 adaptive techniques [Spillane doi:10.1137/15M1028534], or exploitation of the sparsity [Molina and Roux. doi:10.1002/nme.6024]
331
332 Contributed by: P. Gosselet and N. Tardieu
333
334 Example usage:
335 [*] KSP ex71, basic ASM:
336 ./ex71 -pde_type Elasticity -cells 10,10,10 -dim 3 -dirichlet -pc_type asm -mat_type seqaij
337 -ksp_type mpcg -ksp_norm_type unpreconditioned -pc_asm_type basic -pc_asm_overlap 1 -pc_asm_blocks 8
338 -sub_ksp_type preonly -sub_pc_type lu -sub_pc_factor_mat_solver_type mumps -ksp_monitor -ksp_view
339
340 Reference: [Bridson and Greif. doi:10.1137/040620047, Gosselet et al. doi:10.1002/nme.4946]
341
342.seealso: KSPCreate(), KSPSetType(), KSPMPOMIN
343M*/
344PETSC_EXTERNextern __attribute__((visibility ("default"))) PetscErrorCode KSPCreate_MPCG(KSP ksp)
345{
346 PetscErrorCode ierr;
347 KSP_MPCG *mpcg;
348 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
; petscstack->line[petscstack->currentsize] = 348; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);
349 ierr = KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),349,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
350 ierr = KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),350,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
351 ierr = KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),351,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
352 ierr = PetscNewLog(ksp, &mpcg)(PetscMallocA(1,PETSC_TRUE,352,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,(size_t)(1)*sizeof(**(((&mpcg)))),(((&mpcg)))) || PetscLogObjectMemory
((PetscObject)ksp,sizeof(**(&mpcg))));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),352,__func__,"/sandbox/petsc/petsc.next-tmp/src/ksp/ksp/impls/cg/mpcg/mpcg.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
353 ksp->data = (void *)mpcg;
354 ksp->ops->setup = KSPSetUp_MPCG;
355 ksp->ops->solve = KSPSolve_MPCG;
356 ksp->ops->destroy = KSPDestroyDefault;
357 ksp->ops->view = 0;
358 ksp->ops->setfromoptions = KSPSetFromOptions_MPCG;
359 ksp->ops->buildsolution = KSPBuildSolutionDefault;
360 ksp->ops->buildresidual = KSPBuildResidualDefault;
361 PetscFunctionReturn(0)do { do { ; if (petscstack && petscstack->currentsize
> 0) { petscstack->currentsize--; petscstack->function
[petscstack->currentsize] = 0; petscstack->file[petscstack
->currentsize] = 0; petscstack->line[petscstack->currentsize
] = 0; petscstack->petscroutine[petscstack->currentsize
] = PETSC_FALSE; } if (petscstack) { petscstack->hotdepth =
(((petscstack->hotdepth-1)<(0)) ? (0) : (petscstack->
hotdepth-1)); } ; } while (0); return(0);} while (0);
362}