(self, lvl, x, b, cycle)
Multigrid cycling.

        Parameters
        ----------
        lvl : int
            Solve problem on level `lvl`
        x : numpy array
            Initial guess `x` and return correction
        b : numpy array
            Right-hand side for Ax=b
        cycle : {'V','W','F','AMLI'}
            Recursively called cycling function.  The
            Defines the cycling used:
            cycle = 'V',    V-cycle
            cycle = 'W',    W-cycle
            cycle = 'F',    F-cycle
            cycle = 'AMLI', AMLI-cycle
Multigrid cycling.
Multigrid cycling .
def __solve(self, lvl, x, b, cycle):
        """Multigrid cycling.

        Parameters
        ----------
        lvl : int
            Solve problem on level `lvl`
        x : numpy array
            Initial guess `x` and return correction
        b : numpy array
            Right-hand side for Ax=b
        cycle : {'V','W','F','AMLI'}
            Recursively called cycling function.  The
            Defines the cycling used:
            cycle = 'V',    V-cycle
            cycle = 'W',    W-cycle
            cycle = 'F',    F-cycle
            cycle = 'AMLI', AMLI-cycle

        """
        A = self.levels[lvl].A

        self.levels[lvl].presmoother(A, x, b)

        residual = b - A @ x

        coarse_b = self.levels[lvl].R * residual
        coarse_x = np.zeros_like(coarse_b)

        if lvl == len(self.levels) - 2:
            coarse_x[:] = self.coarse_solver(self.levels[-1].A, coarse_b)
        else:
            if cycle == 'V':
                self.__solve(lvl + 1, coarse_x, coarse_b, 'V')
            elif cycle == 'W':
                self.__solve(lvl + 1, coarse_x, coarse_b, cycle)
                self.__solve(lvl + 1, coarse_x, coarse_b, cycle)
            elif cycle == 'F':
                self.__solve(lvl + 1, coarse_x, coarse_b, cycle)
                self.__solve(lvl + 1, coarse_x, coarse_b, 'V')
            elif cycle == 'AMLI':
                # Run nAMLI AMLI cycles, which compute "optimal" corrections by
                # orthogonalizing the coarse-grid corrections in the A-norm
                nAMLI = 2
                Ac = self.levels[lvl + 1].A
                p = np.zeros((nAMLI, coarse_b.shape[0]), dtype=coarse_b.dtype)
                beta = np.zeros((nAMLI, nAMLI), dtype=coarse_b.dtype)
                for k in range(nAMLI):
                    # New search direction --> M^{-1}*residual
                    p[k, :] = 1
                    self.__solve(lvl + 1, p[k, :].reshape(coarse_b.shape),
                                 coarse_b, cycle)

                    # Orthogonalize new search direction to old directions
                    for j in range(k):  # loops from j = 0...(k-1)
                        beta[k, j] = np.inner(p[j, :].conj(), Ac * p[k, :]) /\
                            np.inner(p[j, :].conj(), Ac * p[j, :])
                        p[k, :] -= beta[k, j] * p[j, :]

                    # Compute step size
                    Ap = Ac * p[k, :]
                    alpha = np.inner(p[k, :].conj(), np.ravel(coarse_b)) /\
                        np.inner(p[k, :].conj(), Ap)

                    # Update solution
                    coarse_x += alpha * p[k, :].reshape(coarse_x.shape)

                    # Update residual
                    coarse_b -= alpha * Ap.reshape(coarse_b.shape)
            else:
                raise TypeError(f'Unrecognized cycle type ({cycle})')

        x += self.levels[lvl].P * coarse_x   # coarse grid correction

        self.levels[lvl].postsmoother(A, x, b)
def __solve ( self , lvl , x , b , cycle ) : A = self . levels [ lvl ] . A self . levels [ lvl ] . presmoother ( A , x , b ) residual = b - A @ x coarse_b = self . levels [ lvl ] . R * residual coarse_x = np . zeros_like ( coarse_b ) if lvl == len ( self . levels ) - 2 : coarse_x [ : ] = self . coarse_solver ( self . levels [ - 1 ] . A , coarse_b ) else : if cycle == 'V' : self . __solve ( lvl + 1 , coarse_x , coarse_b , 'V' ) elif cycle == 'W' : self . __solve ( lvl + 1 , coarse_x , coarse_b , cycle ) self . __solve ( lvl + 1 , coarse_x , coarse_b , cycle ) elif cycle == 'F' : self . __solve ( lvl + 1 , coarse_x , coarse_b , cycle ) self . __solve ( lvl + 1 , coarse_x , coarse_b , 'V' ) elif cycle == 'AMLI' : # Run nAMLI AMLI cycles, which compute "optimal" corrections by # orthogonalizing the coarse-grid corrections in the A-norm nAMLI = 2 Ac = self . levels [ lvl + 1 ] . A p = np . zeros ( ( nAMLI , coarse_b . shape [ 0 ] ) , dtype = coarse_b . dtype ) beta = np . zeros ( ( nAMLI , nAMLI ) , dtype = coarse_b . dtype ) for k in range ( nAMLI ) : # New search direction --> M^{-1}*residual p [ k , : ] = 1 self . __solve ( lvl + 1 , p [ k , : ] . reshape ( coarse_b . shape ) , coarse_b , cycle ) # Orthogonalize new search direction to old directions for j in range ( k ) : # loops from j = 0...(k-1) beta [ k , j ] = np . inner ( p [ j , : ] . conj ( ) , Ac * p [ k , : ] ) / np . inner ( p [ j , : ] . conj ( ) , Ac * p [ j , : ] ) p [ k , : ] -= beta [ k , j ] * p [ j , : ] # Compute step size Ap = Ac * p [ k , : ] alpha = np . inner ( p [ k , : ] . conj ( ) , np . ravel ( coarse_b ) ) / np . inner ( p [ k , : ] . conj ( ) , Ap ) # Update solution coarse_x += alpha * p [ k , : ] . reshape ( coarse_x . shape ) # Update residual coarse_b -= alpha * Ap . reshape ( coarse_b . shape ) else : raise TypeError ( f'Unrecognized cycle type ({cycle})' ) x += self . levels [ lvl ] . P * coarse_x # coarse grid correction self . levels [ lvl ] . postsmoother ( A , x , b )