[Scipy-svn] r3844 - in trunk/scipy/sandbox/multigrid: . examples

scipy-svn@scip... scipy-svn@scip...
Wed Jan 16 16:42:53 CST 2008


Author: wnbell
Date: 2008-01-16 16:42:48 -0600 (Wed, 16 Jan 2008)
New Revision: 3844

Modified:
   trunk/scipy/sandbox/multigrid/adaptive.py
   trunk/scipy/sandbox/multigrid/examples/adaptive.py
   trunk/scipy/sandbox/multigrid/multilevel.py
   trunk/scipy/sandbox/multigrid/utils.py
Log:
reworked aSA to better agree with SA


Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2008-01-16 21:57:52 UTC (rev 3843)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2008-01-16 22:42:48 UTC (rev 3844)
@@ -48,148 +48,187 @@
     return As,Ps,Ts,Bs
 
 
-class adaptive_sa_solver:
-    def __init__(self, A, aggregation=None, max_levels=10, max_coarse=100,\
-                 max_candidates=1, mu=5, epsilon=0.1):
+def adaptive_sa_solver(A, max_candidates=1, mu=5, max_levels=10, max_coarse=100, \
+        epsilon=0.0, omega=4.0/3.0, symmetric=True, rescale=True,
+        aggregation=None):
+    """Create a multilevel solver using Smoothed Aggregation (SA)
+
+    *Parameters*:
+
+        A : {csr_matrix}
+            NxN matrix in CSR or BSR format
+        max_candidates : {integer} : default 1
+            Maximum number of near-nullspace candidates to generate.
+        mu : {integer} : default 5
+            Number of cycles used at each level of the adaptive setup phase.
+        max_levels: {integer} : default 10
+            Maximum number of levels to be used in the multilevel solver.
+        max_coarse: {integer} : default 500
+            Maximum number of variables permitted on the coarse grid.
+        epsilon: {float} : default 0.0
+            Strength of connection parameter used in aggregation.
+        omega: {float} : default 4.0/3.0
+            Damping parameter used in prolongator smoothing (0 < omega < 2)
+        symmetric: {boolean} : default True
+            True if A is symmetric, False otherwise
+        rescale: {boolean} : default True
+            If True, symmetrically rescale A by the diagonal
+            i.e. A -> D * A * D,  where D is diag(A)^-0.5
+        aggregation: {None, list of csr_matrix} : optional
+            List of csr_matrix objects that describe a user-defined
+            multilevel aggregation of the variables.
+            TODO ELABORATE
+
+    *Notes*:
+        Unlike the standard Smoothed Aggregation (SA) method, adaptive SA
+        does not require knowledge of near-nullspace candidate vectors.
+        Instead, an adaptive procedure computes one or more candidates 
+        'from scratch'.  This approach is useful when no candidates are known
+        or the candidates have been invalidated due to changes to matrix A.
         
-        self.A  = A
-        self.Rs = []
 
-        if self.A.shape[0] <= max_coarse:
-            raise ValueError,'small matrices not handled yet'
+    *Example*:
+        TODO
 
-        #first candidate
-        x,AggOps = self.__initialization_stage(A, max_levels = max_levels, \
-                                               max_coarse = max_coarse, \
-                                               mu = mu, epsilon = epsilon, \
-                                               aggregation = aggregation )
+    *References*:
+        "Adaptive Smoothed Aggregation ($\alpha$SA) Multigrid"
+            Brezina, Falgout, MacLachlan, Manteuffel, McCormick, and Ruge
+            SIAM Review Volume 47 ,  Issue 2  (2005)
+            http://www.cs.umn.edu/~maclach/research/aSA2.pdf
 
-        #create SA using x here
-        As,Ps,Ts,Bs = sa_hierarchy(A,x,AggOps)
+    """
+    
+    if A.shape[0] <= max_coarse:
+        raise ValueError,'small matrices not handled yet'
 
-        for i in range(max_candidates - 1):
-            x = self.__develop_new_candidate(As,Ps,Ts,Bs,AggOps,mu=mu)
+    #first candidate
+    x,AggOps = asa_initial_setup_stage(A, max_levels = max_levels, \
+            max_coarse = max_coarse, mu = mu, epsilon = epsilon, \
+            aggregation = aggregation )
 
-            B = hstack((Bs[0],x))
-            As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
+    #create SA using x here
+    As,Ps,Ts,Bs = sa_hierarchy(A,x,AggOps)
 
-        #improve candidates?
-        if False:
-            print "improving candidates"
-            B = Bs[0]
-            for i in range(max_candidates):
-                B = B[:,1:]
-                As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
-                x = self.__develop_new_candidate(As,Ps,Ts,Bs,AggOps,mu=mu)
-                B = hstack((B,x))
+    for i in range(max_candidates - 1):
+        x = asa_general_setup_stage(As,Ps,Ts,Bs,AggOps,mu=mu)
+
+        B = hstack((Bs[0],x))
+        As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
+
+    #improve candidates?
+    if False:
+        print "improving candidates"
+        B = Bs[0]
+        for i in range(max_candidates):
+            B = B[:,1:]
             As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
+            x = asa_general_setup_stage(As,Ps,Ts,Bs,AggOps,mu=mu)
+            B = hstack((B,x))
+        As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
 
-        self.Ts = Ts
-        self.solver = multilevel_solver(As,Ps)
-        self.AggOps = AggOps
-        self.Bs = Bs
+    return multilevel_solver(As,Ps)
 
-    def __initialization_stage(self,A,max_levels,max_coarse,mu,epsilon,aggregation):
-        if aggregation is not None:
-            max_coarse = 0
-            max_levels = len(aggregation) + 1
+def asa_initial_setup_stage(A, max_levels, max_coarse, mu, epsilon, aggregation):
+    if aggregation is not None:
+        max_coarse = 0
+        max_levels = len(aggregation) + 1
 
-        # aSA parameters
-        # mu      - number of test relaxation iterations
-        # epsilon - minimum acceptable relaxation convergence factor
+    # aSA parameters
+    # mu      - number of test relaxation iterations
+    # epsilon - minimum acceptable relaxation convergence factor
 
-        #step 1
-        A_l = A
-        x   = rand(A_l.shape[0],1) # TODO see why randn() fails here
-        skip_f_to_i = False
+    #step 1
+    A_l = A
+    x   = rand(A_l.shape[0],1) # TODO see why randn() fails here
+    skip_f_to_i = False
 
-        #step 2
-        gauss_seidel(A_l, x, zeros_like(x), iterations=mu, sweep='symmetric')
+    #step 2
+    gauss_seidel(A_l, x, zeros_like(x), iterations=mu, sweep='symmetric')
 
-        #step 3
-        #TODO test convergence rate here
+    #step 3
+    #TODO test convergence rate here
 
-        As     = [A]
-        Ps     = []
-        AggOps = []
+    As     = [A]
+    Ps     = []
+    AggOps = []
 
-        while len(AggOps) + 1 < max_levels and A_l.shape[0] > max_coarse:
-            if aggregation is None:
-                W_l = sa_constant_interpolation(A_l,epsilon=0) #step 4b
-            else:
-                W_l = aggregation[len(AggOps)]
-            P_l,x = sa_fit_candidates(W_l,x)                   #step 4c
-            I_l   = sa_smoothed_prolongator(A_l,P_l)           #step 4d
-            A_l   = I_l.T.tocsr() * A_l * I_l                  #step 4e
+    while len(AggOps) + 1 < max_levels and A_l.shape[0] > max_coarse:
+        if aggregation is None:
+            W_l = sa_constant_interpolation(A_l,epsilon=0) #step 4b
+        else:
+            W_l = aggregation[len(AggOps)]
+        P_l,x = sa_fit_candidates(W_l,x)                   #step 4c
+        I_l   = sa_smoothed_prolongator(A_l,P_l)           #step 4d
+        A_l   = I_l.T.tocsr() * A_l * I_l                  #step 4e
 
-            AggOps.append(W_l)
-            Ps.append(I_l)
-            As.append(A_l)
+        AggOps.append(W_l)
+        Ps.append(I_l)
+        As.append(A_l)
 
-            if A_l.shape <= max_coarse:  break
+        if A_l.shape <= max_coarse:  break
 
-            if not skip_f_to_i:
-                x_hat = x.copy()                                                   #step 4g
-                gauss_seidel(A_l,x,zeros_like(x),iterations=mu,sweep='symmetric')  #step 4h
-                x_A_x = dot(x.T,A_l*x)
-                if (x_A_x/dot(x_hat.T,A_l*x_hat))**(1.0/mu) < epsilon:             #step 4i
-                    print "sufficient convergence, skipping"
-                    skip_f_to_i = True
-                    if x_A_x == 0:
-                        x = x_hat  #need to restore x
+        if not skip_f_to_i:
+            x_hat = x.copy()                                                   #step 4g
+            gauss_seidel(A_l,x,zeros_like(x),iterations=mu,sweep='symmetric')  #step 4h
+            x_A_x = dot(x.T,A_l*x)
+            if (x_A_x/dot(x_hat.T,A_l*x_hat))**(1.0/mu) < epsilon:             #step 4i
+                print "sufficient convergence, skipping"
+                skip_f_to_i = True
+                if x_A_x == 0:
+                    x = x_hat  #need to restore x
 
-        #update fine-level candidate
-        for A_l,I in reversed(zip(As[1:],Ps)):
-            gauss_seidel(A_l,x,zeros_like(x),iterations=mu,sweep='symmetric')         #TEST
-            x = I * x
-        gauss_seidel(A,x,zeros_like(x),iterations=mu,sweep='symmetric')         #TEST
+    #update fine-level candidate
+    for A_l,I in reversed(zip(As[1:],Ps)):
+        gauss_seidel(A_l,x,zeros_like(x),iterations=mu,sweep='symmetric')         #TEST
+        x = I * x
+    gauss_seidel(A,x,zeros_like(x),iterations=mu,sweep='symmetric')         #TEST
 
-        return x,AggOps  #first candidate,aggregation
+    return x,AggOps  #first candidate,aggregation
 
 
-    def __develop_new_candidate(self,As,Ps,Ts,Bs,AggOps,mu):
-        A = As[0]
+def asa_general_setup_stage(As, Ps, Ts, Bs, AggOps, mu):
+    A = As[0]
 
-        x = randn(A.shape[0],1)
-        b = zeros_like(x)
+    x = randn(A.shape[0],1)
+    b = zeros_like(x)
 
-        x = multilevel_solver(As,Ps).solve(b, x0=x, tol=1e-10, maxiter=mu)
+    x = multilevel_solver(As,Ps).solve(b, x0=x, tol=1e-10, maxiter=mu)
 
-        #TEST FOR CONVERGENCE HERE
+    #TEST FOR CONVERGENCE HERE
 
-        temp_Ps = []
-        temp_As = [A]
+    temp_Ps = []
+    temp_As = [A]
 
-        def make_bridge(P):
-            M,N  = P.shape
-            K    = P.blocksize[0]
-            bnnz = P.indptr[-1]
-            data = zeros( (bnnz, K+1, K), dtype=P.dtype )
-            data[:,:-1,:-1] = P.data
-            return bsr_matrix( (data, P.indices, P.indptr), shape=( (K+1)*(M/K), N) )
+    def make_bridge(P):
+        M,N  = P.shape
+        K    = P.blocksize[0]
+        bnnz = P.indptr[-1]
+        data = zeros( (bnnz, K+1, K), dtype=P.dtype )
+        data[:,:-1,:-1] = P.data
+        return bsr_matrix( (data, P.indices, P.indptr), shape=( (K+1)*(M/K), N) )
 
-        for i in range(len(As) - 2):
-            B = zeros( (x.shape[0], Bs[i+1].shape[1] + 1), dtype=x.dtype)
-            T,R = sa_fit_candidates(AggOp,B)
+    for i in range(len(As) - 2):
+        B = zeros( (x.shape[0], Bs[i+1].shape[1] + 1), dtype=x.dtype)
+        T,R = sa_fit_candidates(AggOp,B)
 
-            P = sa_smoothed_prolongator(A,T)
-            A = P.T.tocsr() * A * P
+        P = sa_smoothed_prolongator(A,T)
+        A = P.T.tocsr() * A * P
 
-            temp_Ps.append(P)
-            temp_As.append(A)
+        temp_Ps.append(P)
+        temp_As.append(A)
 
-            bridge = make_bridge(Ps[i+1])
+        bridge = make_bridge(Ps[i+1])
 
-            solver = multilevel_solver( [A] + As[i+2:], [bridge] + Ps[i+2:] )
+        solver = multilevel_solver( [A] + As[i+2:], [bridge] + Ps[i+2:] )
 
-            x = R[:,-1].reshape(-1,1)
-            x = solver.solve(zeros_like(x), x0=x, tol=1e-8, maxiter=mu)
+        x = R[:,-1].reshape(-1,1)
+        x = solver.solve(zeros_like(x), x0=x, tol=1e-8, maxiter=mu)
 
-        for A,P in reversed(zip(temp_As,temp_Ps)):
-            x = P * x
-            gauss_seidel(A,x,zeros_like(x),iterations=mu,sweep='symmetric')
+    for A,P in reversed(zip(temp_As,temp_Ps)):
+        x = P * x
+        gauss_seidel(A,x,zeros_like(x),iterations=mu,sweep='symmetric')
 
-        return x
+    return x
 
 #    def __augment_cycle(self,As,Ps,Ts,Bs,AggOps,x):
 #        A = As[0]

Modified: trunk/scipy/sandbox/multigrid/examples/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/examples/adaptive.py	2008-01-16 21:57:52 UTC (rev 3843)
+++ trunk/scipy/sandbox/multigrid/examples/adaptive.py	2008-01-16 22:42:48 UTC (rev 3844)
@@ -10,7 +10,7 @@
 
 #A = poisson( (200,200), spacing=(1,1e-2) )  #anisotropic 
 D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
-#A = D * A * D
+A = D * A * D
 
 
 
@@ -28,12 +28,12 @@
 
 print "solving"
 if True:
-    x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=20,tol=1e-12,return_residuals=True)
+    x_sol,residuals = asa.solve(b,x0=x,maxiter=20,tol=1e-12,return_residuals=True)
 else:
     residuals = []
     def add_resid(x):
         residuals.append(linalg.norm(b - A*x))
-    A.psolve = asa.solver.psolve
+    A.psolve = asa.psolve
     x_sol = linalg.cg(A,b,x0=x,maxiter=30,tol=1e-12,callback=add_resid)[0]
 
 residuals = array(residuals)/residuals[0]
@@ -43,7 +43,7 @@
 print "last convergence factor",residuals[-1]/residuals[-2]
 
 print
-print asa.solver
+print asa
 
 print "constant Rayleigh quotient",dot(ones(A.shape[0]),A*ones(A.shape[0]))/float(A.shape[0])
 
@@ -69,10 +69,10 @@
     show()
 
 
-for c in asa.Bs[0].T:
-    plot2d(c)
-    #plot2d_arrows(c)
-    print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
+#for c in asa.Bs[0].T:
+#    plot2d(c)
+#    #plot2d_arrows(c)
+#    print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
 
 #plot2d(x_sol)
 

Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2008-01-16 21:57:52 UTC (rev 3843)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2008-01-16 22:42:48 UTC (rev 3844)
@@ -41,11 +41,9 @@
 
 
 
-def smoothed_aggregation_solver(A, B=None,  \
-                                aggregation=None, max_levels=10, \
-                                max_coarse=500, epsilon=0.0, \
-                                omega=4.0/3.0, symmetric=True, \
-                                rescale = True):
+def smoothed_aggregation_solver(A, B=None, max_levels=10, max_coarse=500, \
+        epsilon=0.0, omega=4.0/3.0, symmetric=True, rescale=True, \
+        aggregation=None):
     """Create a multilevel solver using Smoothed Aggregation (SA)
 
     *Parameters*:
@@ -55,19 +53,6 @@
         B : {None, array_like} : optional
             Near-nullspace candidates stored in the columns of an NxK array.
             The default value B=None is equivalent to B=ones((N,1))
-        blocks : {None, array_like} : optional
-            Array of length N that groups the variables into 'superblocks'.
-            For example, in a 2d vector-valued problem where the even
-            variables [0,2,4,...N-2] correspond to the x-components and the
-            odd variables [1,3,5,...,N-1] correspond to the y-components then
-            blocks=[0,0,1,1,2,2,...,N/2,N/2] is expected.  The default
-            value blocks=None is equivalent to blocks=[0,1,2,..,N] which
-            implies that each variable should be aggregated seperately.
-            The default is appropriate for scalar valued problems.
-        aggregation: {None, list of csr_matrix} : optional
-            List of csr_matrix objects that describe a user-defined
-            multilevel aggregation of the variables.
-            TODO ELABORATE
         max_levels: {integer} : default 10
             Maximum number of levels to be used in the multilevel solver.
         max_coarse: {integer} : default 500
@@ -81,6 +66,10 @@
         rescale: {boolean} : default True
             If True, symmetrically rescale A by the diagonal
             i.e. A -> D * A * D,  where D is diag(A)^-0.5
+        aggregation: {None, list of csr_matrix} : optional
+            List of csr_matrix objects that describe a user-defined
+            multilevel aggregation of the variables.
+            TODO ELABORATE
 
     *Example*:
         TODO

Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py	2008-01-16 21:57:52 UTC (rev 3843)
+++ trunk/scipy/sandbox/multigrid/utils.py	2008-01-16 22:42:48 UTC (rev 3844)
@@ -21,6 +21,7 @@
 
         tol : {scalar}
             Tolerance of approximation
+            Currently unused
 
         maxiter : {integer}
             Maximum number of iterations to perform



More information about the Scipy-svn mailing list