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

scipy-svn@scip... scipy-svn@scip...
Tue Jan 15 07:54:57 CST 2008


Author: wnbell
Date: 2008-01-15 07:54:53 -0600 (Tue, 15 Jan 2008)
New Revision: 3837

Added:
   trunk/scipy/sandbox/multigrid/examples/
   trunk/scipy/sandbox/multigrid/examples/adaptive.py
   trunk/scipy/sandbox/multigrid/examples/dec_example.py
   trunk/scipy/sandbox/multigrid/examples/simple_example.py
Modified:
   trunk/scipy/sandbox/multigrid/adaptive.py
Log:
updating adaptive SA code to use BSR


Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2008-01-15 13:09:03 UTC (rev 3836)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2008-01-15 13:54:53 UTC (rev 3837)
@@ -1,4 +1,5 @@
 import numpy,scipy,scipy.sparse
+
 from numpy import sqrt, ravel, diff, zeros, zeros_like, inner, concatenate, \
                   asarray, hstack, ascontiguousarray, isinf, dot
 from numpy.random import randn
@@ -11,89 +12,8 @@
 
 
 
-def augment_candidates(AggOp, old_Q, old_R, new_candidate):
-    #TODO update P and A also
 
-    K = old_R.shape[1]
 
-    #determine blocksizes
-    if new_candidate.shape[0] == old_Q.shape[0]:
-        #then this is the first prolongator
-        old_bs = (1,K)
-        new_bs = (1,K+1)
-    else:
-        old_bs = (K,K)
-        new_bs = (K+1,K+1)
-
-    AggOp = expand_into_blocks(AggOp,new_bs[0],1).tocsr() #TODO switch to block matrix
-
-
-    # tentative prolongator
-    #TODO USE BSR
-    Q_indptr  = (K+1)*AggOp.indptr
-    Q_indices = ((K+1)*AggOp.indices).repeat(K+1)
-    for i in range(K+1):
-        Q_indices[i::K+1] += i
-    Q_data = zeros((AggOp.indptr[-1]/new_bs[0],) + new_bs)
-    Q_data[:,:old_bs[0],:old_bs[1]] = old_Q.data.reshape((-1,) + old_bs)   #TODO BSR change
-
-    # coarse candidates
-    R = zeros((AggOp.shape[1],K+1,K+1))
-    R[:,:K,:K] = old_R.reshape(-1,K,K)
-
-    c = new_candidate.reshape(-1)[diff(AggOp.indptr) == 1]  #eliminate DOFs that aggregation misses
-    threshold = 1e-10 * abs(c).max()   # cutoff for small basis functions
-
-    X = csr_matrix((c,AggOp.indices,AggOp.indptr),shape=AggOp.shape)
-
-    #orthogonalize X against previous
-    for i in range(K):
-        old_c = ascontiguousarray(Q_data[:,:,i].reshape(-1))
-        D_AtX = csr_matrix((old_c*X.data,X.indices,X.indptr),shape=X.shape).sum(axis=0).A.flatten() #same as diagonal of A.T * X
-        R[:,i,K] = D_AtX
-        X.data -= D_AtX[X.indices] * old_c
-
-    #normalize X
-    D_XtX = csr_matrix((X.data**2,X.indices,X.indptr),shape=X.shape).sum(axis=0).A.flatten() #same as diagonal of X.T * X
-    col_norms = sqrt(D_XtX)
-    mask = col_norms < threshold  # find small basis functions
-    col_norms[mask] = 0           # and set them to zero
-
-    R[:,K,K] = col_norms      # store diagonal entry into R
-
-    col_norms = 1.0/col_norms
-    col_norms[mask] = 0
-    X.data *= col_norms[X.indices]
-    Q_data[:,:,-1] = X.data.reshape(-1,new_bs[0])
-
-    Q_data = Q_data.reshape(-1)  #TODO BSR change
-    R = R.reshape(-1,K+1)
-
-    Q = csr_matrix((Q_data,Q_indices,Q_indptr),shape=(AggOp.shape[0],(K+1)*AggOp.shape[1]))
-
-    return Q,R
-
-
-
-
-
-
-def smoothed_prolongator(P,A):
-    #just use Richardson for now
-    #omega = 4.0/(3.0*approximate_spectral_radius(A))
-    #return P - omega*(A*P)
-    #return P  #TEST
-
-    D = diag_sparse(A)
-    D_inv_A = diag_sparse(1.0/D)*A
-    omega = 4.0/(3.0*approximate_spectral_radius(D_inv_A))
-    print "spectral radius",approximate_spectral_radius(D_inv_A) #TODO remove this
-    D_inv_A *= omega
-
-    return P - D_inv_A*P
-
-
-
 def sa_hierarchy(A,B,AggOps):
     """
     Construct multilevel hierarchy using Smoothed Aggregation
@@ -114,7 +34,7 @@
 
     for AggOp in AggOps:
         P,B = sa_fit_candidates(AggOp,B)
-        I   = smoothed_prolongator(P,A)
+        I   = sa_smoothed_prolongator(P,A)
         A   = I.T.tocsr() * A * I
         As.append(A)
         Ts.append(P)
@@ -148,13 +68,13 @@
             x = self.__develop_new_candidate(As,Ps,Ts,Bs,AggOps,mu=mu)
 
             #TODO which is faster?
-            As,Ps,Ts,Bs = self.__augment_cycle(As,Ps,Ts,Bs,AggOps,x)
+            #As,Ps,Ts,Bs = self.__augment_cycle(As,Ps,Ts,Bs,AggOps,x)
 
-            #B = hstack((Bs[0],x))
-            #As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
+            B = hstack((Bs[0],x))
+            As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
 
         #improve candidates?
-        if True:
+        if False:
             print "improving candidates"
             B = Bs[0]
             for i in range(max_candidates):
@@ -199,7 +119,7 @@
             else:
                 W_l   = aggregation[len(AggOps)]
             P_l,x = sa_fit_candidates(W_l,x)                             #step 4c
-            I_l   = smoothed_prolongator(P_l,A_l)                          #step 4d
+            I_l   = sa_smoothed_prolongator(P_l,A_l)                          #step 4d
             A_l   = I_l.T.tocsr() * A_l * I_l                              #step 4e
 
             blocks = None #not needed on subsequent levels
@@ -244,24 +164,25 @@
         temp_Ps = []
         temp_As = [A]
 
-        def make_bridge(P,K):
-            indptr = P.indptr[:-1].reshape(-1,K-1)
-            indptr = hstack((indptr,indptr[:,-1].reshape(-1,1)))
-            indptr = indptr.reshape(-1)
-            indptr = hstack((indptr,indptr[-1:])) #duplicate last element
-            return csr_matrix((P.data,P.indices,indptr),shape=(K*P.shape[0]/(K-1),P.shape[1]))
+        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):
-            T,R = augment_candidates(AggOps[i], Ts[i], Bs[i+1], x)
+            B = zeros( (x.shape[0], Bs[i+1].shape[1] + 1), dtype=x.dtype)
+            T,R = sa_fit_candidates(AggOp,B)
 
-            P = smoothed_prolongator(T,A)
+            P = sa_smoothed_prolongator(T,A)
             A = P.T.tocsr() * A * P
 
             temp_Ps.append(P)
             temp_As.append(A)
 
-            #TODO USE BSR (K,K) -> (K,K-1)
-            bridge = make_bridge(Ps[i+1],R.shape[1])
+            bridge = make_bridge(Ps[i+1])
 
             solver = multilevel_solver( [A] + As[i+2:], [bridge] + Ps[i+2:] )
 
@@ -274,111 +195,87 @@
 
         return x
 
-    def __augment_cycle(self,As,Ps,Ts,Bs,AggOps,x):
-        A = As[0]
+#    def __augment_cycle(self,As,Ps,Ts,Bs,AggOps,x):
+#        A = As[0]
+#
+#        new_As = [A]
+#        new_Ps = []
+#        new_Ts = []
+#        new_Bs = [ hstack((Bs[0],x)) ]
+#
+#        for i in range(len(As) - 1):
+#            T,R = augment_candidates(AggOps[i], Ts[i], Bs[i+1], x)
+#
+#            P = sa_smoothed_prolongator(T,A)
+#            A = P.T.tocsr() * A * P
+#
+#            new_As.append(A)
+#            new_Ps.append(P)
+#            new_Ts.append(T)
+#            new_Bs.append(R)
+#
+#            x = R[:,-1].reshape(-1,1)
+#
+#        return new_As,new_Ps,new_Ts,new_Bs
 
-        new_As = [A]
-        new_Ps = []
-        new_Ts = []
-        new_Bs = [ hstack((Bs[0],x)) ]
 
-        for i in range(len(As) - 1):
-            T,R = augment_candidates(AggOps[i], Ts[i], Bs[i+1], x)
 
-            P = smoothed_prolongator(T,A)
-            A = P.T.tocsr() * A * P
 
-            new_As.append(A)
-            new_Ps.append(P)
-            new_Ts.append(T)
-            new_Bs.append(R)
 
-            x = R[:,-1].reshape(-1,1)
+#def augment_candidates(AggOp, old_Q, old_R, new_candidate):
+#    #TODO update P and A also
+#
+#    K = old_R.shape[1]
+#
+#    #determine blocksizes
+#    if new_candidate.shape[0] == old_Q.shape[0]:
+#        #then this is the first prolongator
+#        old_bs = (1,K)
+#        new_bs = (1,K+1)
+#    else:
+#        old_bs = (K,K)
+#        new_bs = (K+1,K+1)
+#
+#    #AggOp = expand_into_blocks(AggOp,new_bs[0],1).tocsr() #TODO switch to block matrix
+#
+#    # tentative prolongator
+#    Q_data = zeros( (AggOp.indptr[-1],) + new_bs)
+#    Q_data[:,:old_bs[0],:old_bs[1]] = old_Q.data.reshape((-1,) + old_bs)
+#
+#    # coarse candidates
+#    R = zeros((AggOp.shape[1],K+1,K+1))
+#    R[:,:K,:K] = old_R.reshape(-1,K,K)
+#
+#    c = new_candidate.reshape(-1)[diff(AggOp.indptr) == 1]  #eliminate DOFs that aggregation misses
+#    threshold = 1e-10 * abs(c).max()   # cutoff for small basis functions
+#
+#    X = bsr_matrix((c,AggOp.indices,AggOp.indptr),shape=AggOp.shape)
+#
+#    #orthogonalize X against previous
+#    for i in range(K):
+#        old_c = ascontiguousarray(Q_data[:,:,i].reshape(-1))
+#        D_AtX = csr_matrix((old_c*X.data,X.indices,X.indptr),shape=X.shape).sum(axis=0).A.flatten() #same as diagonal of A.T * X
+#        R[:,i,K] = D_AtX
+#        X.data -= D_AtX[X.indices] * old_c
+#
+#    #normalize X
+#    D_XtX = csr_matrix((X.data**2,X.indices,X.indptr),shape=X.shape).sum(axis=0).A.flatten() #same as diagonal of X.T * X
+#    col_norms = sqrt(D_XtX)
+#    mask = col_norms < threshold  # find small basis functions
+#    col_norms[mask] = 0           # and set them to zero
+#
+#    R[:,K,K] = col_norms      # store diagonal entry into R
+#
+#    col_norms = 1.0/col_norms
+#    col_norms[mask] = 0
+#    X.data *= col_norms[X.indices]
+#    Q_data[:,:,-1] = X.data.reshape(-1,new_bs[0])
+#
+#    Q_data = Q_data.reshape(-1)  #TODO BSR change
+#    R = R.reshape(-1,K+1)
+#
+#    Q = csr_matrix((Q_data,Q_indices,Q_indptr),shape=(AggOp.shape[0],(K+1)*AggOp.shape[1]))
+#
+#    return Q,R
 
-        return new_As,new_Ps,new_Ts,new_Bs
 
-
-if __name__ == '__main__':
-    from scipy import *
-    from utils import diag_sparse
-    from multilevel import poisson_problem1D,poisson_problem2D
-
-    blocks = None
-    aggregation = None
-
-    #A = poisson_problem2D(200,1e-2)
-    #aggregation = [ sa_constant_interpolation(A*A*A,epsilon=0.0) ]
-
-    #A = io.mmread("tests/sample_data/laplacian_41_3dcube.mtx").tocsr()
-    #A = io.mmread("laplacian_40_3dcube.mtx").tocsr()
-    #A = io.mmread("/home/nathan/Desktop/9pt/9pt-100x100.mtx").tocsr()
-    #A = io.mmread("/home/nathan/Desktop/BasisShift_W_EnergyMin_Luke/9pt-5x5.mtx").tocsr()
-
-
-    #D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
-    #A = D * A * D
-
-    A = io.mmread("tests/sample_data/elas30_A.mtx").tocsr()
-    blocks = arange(A.shape[0]/2).repeat(2)
-
-    from time import clock; start = clock()
-    asa = adaptive_sa_solver(A,max_candidates=3,mu=5,blocks=blocks,aggregation=aggregation)
-    print "Adaptive Solver Construction: %s seconds" % (clock() - start); del start
-
-    scipy.random.seed(0)  #make tests repeatable
-    x = randn(A.shape[0])
-    b = A*randn(A.shape[0])
-    #b = zeros(A.shape[0])
-
-
-    print "solving"
-    if False:
-        x_sol,residuals = asa.solver.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
-        x_sol = linalg.cg(A,b,x0=x,maxiter=30,tol=1e-12,callback=add_resid)[0]
-
-    residuals = array(residuals)/residuals[0]
-
-    print "residuals ",residuals
-    print "mean convergence factor",(residuals[-1]/residuals[0])**(1.0/len(residuals))
-    print "last convergence factor",residuals[-1]/residuals[-2]
-
-    print
-    print asa.solver
-
-    print "constant Rayleigh quotient",dot(ones(A.shape[0]),A*ones(A.shape[0]))/float(A.shape[0])
-
-    def plot2d_arrows(x):
-        from pylab import figure,quiver,show
-        x = x.reshape(-1)
-        N = (len(x)/2)**0.5
-        assert(2 * N * N == len(x))
-        X = linspace(-1,1,N).reshape(1,N).repeat(N,0).reshape(-1)
-        Y = linspace(-1,1,N).reshape(1,N).repeat(N,0).T.reshape(-1)
-
-        dX = x[0::2]
-        dY = x[1::2]
-
-        figure()
-        quiver(X,Y,dX,dY)
-        show()
-
-    def plot2d(x):
-        from pylab import pcolor,figure,show
-        figure()
-        pcolor(x.reshape(sqrt(len(x)),sqrt(len(x))))
-        show()
-
-
-    for c in asa.Bs[0].T:
-        #plot2d(c)
-        plot2d_arrows(c)
-        print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
-
-
-    ##W = asa.AggOps[0]*asa.AggOps[1]
-    ##pcolor((W * rand(W.shape[1])).reshape((200,200)))

Added: trunk/scipy/sandbox/multigrid/examples/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/examples/adaptive.py	2008-01-15 13:09:03 UTC (rev 3836)
+++ trunk/scipy/sandbox/multigrid/examples/adaptive.py	2008-01-15 13:54:53 UTC (rev 3837)
@@ -0,0 +1,84 @@
+
+from scipy import *
+from scipy.sandbox.multigrid.utils import diag_sparse
+from scipy.sandbox.multigrid.gallery import poisson, linear_elasticity
+
+
+#A = poisson( (200,200), spacing=(1,1e-2)
+#aggregation = [ sa_constant_interpolation(A*A*A,epsilon=0.0) ]
+
+#A = io.mmread("tests/sample_data/laplacian_41_3dcube.mtx").tocsr()
+#A = io.mmread("laplacian_40_3dcube.mtx").tocsr()
+#A = io.mmread("/home/nathan/Desktop/9pt/9pt-100x100.mtx").tocsr()
+#A = io.mmread("/home/nathan/Desktop/BasisShift_W_EnergyMin_Luke/9pt-5x5.mtx").tocsr()
+
+
+#D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
+#A = D * A * D
+
+#A = io.mmread("tests/sample_data/elas30_A.mtx").tocsr()
+
+A,B = linear_elasticity( (100,100) )
+
+from time import clock; start = clock()
+asa = adaptive_sa_solver(A,max_candidates=3,mu=5,blocks=blocks,aggregation=aggregation)
+print "Adaptive Solver Construction: %s seconds" % (clock() - start); del start
+
+scipy.random.seed(0)  #make tests repeatable
+x = randn(A.shape[0])
+b = A*randn(A.shape[0])
+#b = zeros(A.shape[0])
+
+
+print "solving"
+if False:
+    x_sol,residuals = asa.solver.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
+    x_sol = linalg.cg(A,b,x0=x,maxiter=30,tol=1e-12,callback=add_resid)[0]
+
+residuals = array(residuals)/residuals[0]
+
+print "residuals ",residuals
+print "mean convergence factor",(residuals[-1]/residuals[0])**(1.0/len(residuals))
+print "last convergence factor",residuals[-1]/residuals[-2]
+
+print
+print asa.solver
+
+print "constant Rayleigh quotient",dot(ones(A.shape[0]),A*ones(A.shape[0]))/float(A.shape[0])
+
+def plot2d_arrows(x):
+    from pylab import figure,quiver,show
+    x = x.reshape(-1)
+    N = (len(x)/2)**0.5
+    assert(2 * N * N == len(x))
+    X = linspace(-1,1,N).reshape(1,N).repeat(N,0).reshape(-1)
+    Y = linspace(-1,1,N).reshape(1,N).repeat(N,0).T.reshape(-1)
+
+    dX = x[0::2]
+    dY = x[1::2]
+
+    figure()
+    quiver(X,Y,dX,dY)
+    show()
+
+def plot2d(x):
+    from pylab import pcolor,figure,show
+    figure()
+    pcolor(x.reshape(sqrt(len(x)),sqrt(len(x))))
+    show()
+
+
+for c in asa.Bs[0].T:
+    #plot2d(c)
+    plot2d_arrows(c)
+    print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
+
+
+##W = asa.AggOps[0]*asa.AggOps[1]
+##pcolor((W * rand(W.shape[1])).reshape((200,200)))
+

Added: trunk/scipy/sandbox/multigrid/examples/dec_example.py
===================================================================
--- trunk/scipy/sandbox/multigrid/examples/dec_example.py	2008-01-15 13:09:03 UTC (rev 3836)
+++ trunk/scipy/sandbox/multigrid/examples/dec_example.py	2008-01-15 13:54:53 UTC (rev 3837)
@@ -0,0 +1,241 @@
+
+from scipy import *
+from scipy.sparse import *
+from pydec import *
+from pydec.multigrid.discrete_laplacian import boundary_hierarchy, discrete_laplacian_solver, hodge_solver
+
+from scipy.sandbox.multigrid import smoothed_aggregation_solver,multigridtools,multilevel_solver
+from scipy.sandbox.multigrid.adaptive import adaptive_sa_solver
+from scipy.sandbox.multigrid.sa import sa_smoothed_prolongator
+from scipy.sandbox.multigrid.utils import expand_into_blocks
+
+
+## Load mesh from file
+mesh_path = '../../../../../hirani_group/wnbell/meshes/'
+#mesh = read_mesh(mesh_path + 'rocket/rocket.xml')
+#mesh = read_mesh(mesh_path + 'genus3/genus3_168k.xml')
+#mesh = read_mesh(mesh_path + 'genus3/genus3_455k.xml')
+#mesh = read_mesh(mesh_path + '/torus/torus.xml')
+mesh = read_mesh(mesh_path + '/sq14tri/sq14tri.xml')
+for i in range(5):
+    mesh['vertices'],mesh['elements'] = loop_subdivision(mesh['vertices'],mesh['elements'])
+cmplx = simplicial_complex(mesh['vertices'],mesh['elements'])
+
+## Construct mesh manually
+#bitmap = ones((60,60),dtype='bool')
+#bitmap[1::10,1::10] = False
+#bitmap[100:150,100:400] = False
+#cmplx = regular_cube_complex(regular_cube_mesh(bitmap))
+
+def curl_curl_prolongator(D_nodal,vertices):
+    if not isspmatrix_csr(D_nodal):
+        raise TypeError('expected csr_matrix')
+
+    A = D_nodal.T.tocsr() * D_nodal
+    aggs = multigridtools.sa_get_aggregates(A.shape[0],A.indptr,A.indices)
+
+    num_edges = D_nodal.shape[0]
+    num_basis = vertices.shape[1]
+    num_aggs  = aggs.max() + 1
+
+    # replace with CSR + eliminate duplicates
+    #indptr  = (2*num_basis) * arange(num_edges+1)
+    ## same same
+    #csr_matrix((data,indices,indptr),shape=(num_edges,num_aggs))
+
+    row  = arange(num_edges).repeat(2*num_basis)
+    col  = (num_basis*aggs[D_nodal.indices]).repeat(num_basis)
+    col = col.reshape(-1,num_basis) + arange(num_basis)
+    col = col.reshape(-1)
+    data = tile(0.5 * (D_nodal*vertices),(1,2)).reshape(-1)
+
+    return coo_matrix((data,(row,col)),shape=(num_edges,num_basis*num_aggs)).tocsr()
+
+
+
+
+
+def whitney_innerproduct_cache(cmplx,k):
+    h = hash(cmplx.vertices.tostring()) ^ hash(cmplx.simplices.tostring()) ^ hash(k)
+
+    filename = "/home/nathan/.pydec/cache/whitney_" + str(h) + ".mtx"
+
+    try:
+        import pydec
+        M = pydec.io.read_array(filename)
+    except:
+        import pydec
+        M = whitney_innerproduct(cmplx,k)
+        pydec.io.write_array(filename,M)
+
+    return M
+
+
+
+def cube_innerproduct_cache(cmplx,k):
+    h = hash(cmplx.mesh.bitmap.tostring()) ^ hash(cmplx.mesh.bitmap.shape) ^ hash(k)
+
+    filename = "/home/nathan/.pydec/cache/cube_" + str(h) + ".mtx"
+
+    try:
+        import pydec
+        M = pydec.io.read_array(filename)
+    except:
+        import pydec
+        M = regular_cube_innerproduct(cmplx,k)
+        pydec.io.write_array(filename,M)
+
+    return M
+
+
+
+#solve d_k d_k problem for all reasonable k
+#from pylab import semilogy,show,xlabel,ylabel,legend,ylim,xlim
+#from matplotlib.font_manager import fontManager, FontProperties
+
+cochain_complex = cmplx.cochain_complex()
+
+for i in [1]: #range(len(cochain_complex)-1):
+    print "computing mass matrix"
+
+    if isinstance(cmplx,simplicial_complex):
+        Mi = whitney_innerproduct_cache(cmplx,i+1)
+    else:
+        Mi = regular_cube_innerproduct(cmplx,i+1)
+
+
+    dimension = mesh['vertices'].shape[1]
+
+    if True:
+
+        d0 = cmplx[0].d
+        d1 = cmplx[1].d
+
+        #A = (d1.T.tocsr() * d1 + d0 * d0.T.tocsr()).astype('d')
+        A = (d1.T.tocsr()  * d1).astype('d')
+
+        P = curl_curl_prolongator(d0,mesh['vertices'])
+
+        num_blocks = P.shape[1]/dimension
+        blocks = arange(num_blocks).repeat(dimension)
+
+        P = sa_smoothed_prolongator(A,P,epsilon=0,omega=4.0/3.0)
+
+        PAP = P.T.tocsr() * A * P
+
+        candidates = None
+        candidates = zeros((num_blocks,dimension,dimension))
+        for n in range(dimension):
+            candidates[:,n,n] = 1.0
+        candidates = candidates.reshape(-1,dimension)
+
+        ml = smoothed_aggregation_solver(PAP,epsilon=0.0,candidates=candidates,blocks=blocks)
+        #A = PAP
+        ml = multilevel_solver([A] + ml.As, [P] + ml.Ps)
+    else:
+
+        bh = boundary_hierarchy(cochain_complex)
+        while len(bh) < 3:
+            bh.coarsen()
+        print repr(bh)
+
+        N = len(cochain_complex) - 1
+
+        B =  bh[0][N - i].B
+
+        A = (B.T.tocsr() * B).astype('d')
+        #A = B.T.tocsr() * Mi * B
+
+        constant_prolongators = [lvl[N - i].I for lvl in bh[:-1]]
+
+        method = 'aSA'
+
+        if method == 'RS':
+            As = [A]
+            Ps = []
+            for T in constant_prolongators:
+                Ps.append( sa_smoothed_prolongator(As[-1],T,epsilon=0.0,omega=4.0/3.0) )
+                As.append(Ps[-1].T.tocsr() * As[-1] * Ps[-1])
+            ml = multilevel_solver(As,Ps)
+
+        else:
+            if method == 'BSA':
+                if i == 0:
+                    candidates = None
+                else:
+                    candidates = cmplx[0].d * mesh['vertices']
+                    K = candidates.shape[1]
+
+                    constant_prolongators = [constant_prolongators[0]] + \
+                            [expand_into_blocks(T,K,1).tocsr() for T in constant_prolongators[1:] ]
+
+                    ml = smoothed_aggregation_solver(A,candidates,aggregation=constant_prolongators)
+            elif method == 'aSA':
+                asa = adaptive_sa_solver(A,aggregation=constant_prolongators,max_candidates=dimension,epsilon=0.0)
+                ml = asa.solver
+            else:
+                raise ValuerError,'unknown method'
+
+        #ml = smoothed_aggregation_solver(A,candidates)
+
+    #x = d0 * mesh['vertices'][:,0]
+    x = rand(A.shape[0])
+    b = zeros_like(x)
+    #b = A*rand(A.shape[0])
+
+    if True:
+        x_sol,residuals = ml.solve(b,x0=x,maxiter=50,tol=1e-12,return_residuals=True)
+    else:
+        residuals = []
+        def add_resid(x):
+            residuals.append(linalg.norm(b - A*x))
+        A.psolve = ml.psolve
+
+        from pydec import cg
+        x_sol = cg(A,b,x0=x,maxiter=40,tol=1e-8,callback=add_resid)[0]
+
+
+    residuals = array(residuals)/residuals[0]
+    avg_convergence_ratio = residuals[-1]**(1.0/len(residuals))
+    print "average convergence ratio",avg_convergence_ratio
+    print "last convergence ratio",residuals[-1]/residuals[-2]
+
+    print residuals
+
+
+
+
+##candidates = None
+##blocks = None
+##
+##
+##
+##A = io.mmread('tests/sample_data/elas30_A.mtx').tocsr()
+##candidates = io.mmread('tests/sample_data/elas30_nullspace.mtx')
+##candidates = [ array(candidates[:,x]) for x in range(candidates.shape[1]) ]
+##blocks = arange(A.shape[0]/2).repeat(2)
+##
+##ml = smoothed_aggregation_solver(A,candidates,blocks=blocks,epsilon=0,max_coarse=10,max_levels=10)
+###ml = ruge_stuben_solver(A)
+##
+##x = rand(A.shape[0])
+###b = zeros_like(x)
+##b = A*rand(A.shape[0])
+##
+##if True:
+##    x_sol,residuals = ml.solve(b,x0=x,maxiter=30,tol=1e-12,return_residuals=True)
+##else:
+##    residuals = []
+##    def add_resid(x):
+##        residuals.append(linalg.norm(b - A*x))
+##    A.psolve = ml.psolve
+##    x_sol = linalg.cg(A,b,x0=x,maxiter=25,tol=1e-12,callback=add_resid)[0]
+##
+##
+##residuals = array(residuals)/residuals[0]
+##avg_convergence_ratio = residuals[-1]**(1.0/len(residuals))
+##print "average convergence ratio",avg_convergence_ratio
+##print "last convergence ratio",residuals[-1]/residuals[-2]
+##
+##print residuals
+##

Added: trunk/scipy/sandbox/multigrid/examples/simple_example.py
===================================================================
--- trunk/scipy/sandbox/multigrid/examples/simple_example.py	2008-01-15 13:09:03 UTC (rev 3836)
+++ trunk/scipy/sandbox/multigrid/examples/simple_example.py	2008-01-15 13:54:53 UTC (rev 3837)
@@ -0,0 +1,28 @@
+from scipy import *
+from scipy.sandbox.multigrid.sa import *
+from scipy.sandbox.multigrid import *
+from scipy.sandbox.multigrid.utils import *
+from time import clock
+
+mats = io.loadmat('/home/nathan/Work/ogroup/matrices/elasticity/simple_grid_2d/elasticity_500x500.mat')
+A = mats['A'].tobsr(blocksize=(2,2))
+B = mats['B']
+
+#A = poisson_problem2D(50)
+#B = None
+
+start = clock()
+sa = smoothed_aggregation_solver(A,B=B)
+print "constructed solver in %s seconds" % (clock() - start)
+
+x0 = rand(A.shape[0])
+b = zeros_like(x0)
+start = clock()
+x,residuals = sa.solve(b,x0=x0,return_residuals=True)
+print "solved in %s seconds" % (clock() - start)
+    
+residuals = array(residuals)/residuals[0]
+avg_convergence_ratio = residuals[-1]**(1.0/len(residuals))
+print "average convergence ratio",avg_convergence_ratio
+print "last convergence ratio",residuals[-1]/residuals[-2]
+print 'residuals',residuals



More information about the Scipy-svn mailing list