# [Scipy-svn] r3478 - trunk/scipy/sandbox/multigrid

scipy-svn@scip... scipy-svn@scip...
Wed Oct 31 20:31:40 CDT 2007

```Author: wnbell
Date: 2007-10-31 20:31:37 -0500 (Wed, 31 Oct 2007)
New Revision: 3478

Modified:
trunk/scipy/sandbox/multigrid/dec_test.py
trunk/scipy/sandbox/multigrid/multilevel.py
trunk/scipy/sandbox/multigrid/sa.py
Log:
added new method for curl-curl problem

===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-31 21:46:29 UTC (rev 3477)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2007-11-01 01:31:37 UTC (rev 3478)
@@ -1,18 +1,18 @@
import numpy,scipy,scipy.sparse
from numpy import sqrt, ravel, diff, zeros, zeros_like, inner, concatenate, \
-                  asarray, hstack, ascontiguousarray, isinf
+                  asarray, hstack, ascontiguousarray, isinf, dot
from numpy.random import randn
from scipy.sparse import csr_matrix,coo_matrix

from relaxation import gauss_seidel
from multilevel import multilevel_solver
from sa import sa_constant_interpolation,sa_fit_candidates

def augment_candidates(AggOp, old_Q, old_R, new_candidate):
-    #TODO update P also
+    #TODO update P and A also

K = old_R.shape[1]

@@ -124,7 +124,7 @@

-    def __init__(self, A, blocks=None, max_levels=10, max_coarse=100,\
+    def __init__(self, A, blocks=None, aggregation=None, max_levels=10, max_coarse=100,\
max_candidates=1, mu=5, epsilon=0.1):

self.A = A
@@ -138,8 +138,9 @@
x,AggOps = self.__initialization_stage(A, blocks = blocks, \
max_levels = max_levels, \
max_coarse = max_coarse, \
-                                               mu = mu, epsilon = epsilon)
-
+                                               mu = mu, epsilon = epsilon, \
+                                               aggregation = aggregation )
+
#create SA using x here
As,Ps,Ts,Bs = sa_hierarchy(A,x,AggOps)

@@ -147,8 +148,20 @@
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)
-            B = hstack((Bs[0],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)
+
+        #improve candidates?
+        if True:
+            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))
As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)

self.Ts = Ts
@@ -156,12 +169,12 @@
self.AggOps = AggOps
self.Bs = Bs

+    def __initialization_stage(self,A,blocks,max_levels,max_coarse,mu,epsilon,aggregation):
+        if aggregation is not None:
+            max_coarse = 0
+            max_levels = len(aggregation) + 1

-    def __initialization_stage(self,A,blocks,max_levels,max_coarse,mu,epsilon):
-        AggOps = []
-        Ps     = []
-
-        # aSA parameters
+        # aSA parameters
# mu      - number of test relaxation iterations
# epsilon - minimum acceptable relaxation convergence factor

@@ -177,9 +190,14 @@
#TODO test convergence rate here

As = [A]
+        AggOps = []
+        Ps = []

-        while len(AggOps) + 1 < max_levels  and A_l.shape[0] > max_coarse:
-            W_l   = sa_constant_interpolation(A_l,epsilon=0,blocks=blocks) #step 4b
+        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,blocks=blocks) #step 4b
+            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
A_l   = I_l.T.tocsr() * A_l * I_l                              #step 4e
@@ -286,8 +304,11 @@
from multilevel import poisson_problem1D,poisson_problem2D

blocks = None
+    aggregation = None

-    A = poisson_problem2D(100)
+    #A = poisson_problem2D(200,1e-2)
+    #aggregation = [ sa_constant_interpolation(A*A*A,epsilon=0.0) ]
+
@@ -301,17 +322,17 @@
blocks = arange(A.shape[0]/2).repeat(2)

from time import clock; start = clock()
print "Adaptive Solver Construction: %s seconds" % (clock() - start); del start

-    #scipy.random.seed(0)  #make tests repeatable
+    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 True:
+    if False:
x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=20,tol=1e-12,return_residuals=True)
else:
residuals = []
@@ -352,11 +373,12 @@
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)))

Modified: trunk/scipy/sandbox/multigrid/dec_test.py
===================================================================
--- trunk/scipy/sandbox/multigrid/dec_test.py	2007-10-31 21:46:29 UTC (rev 3477)
+++ trunk/scipy/sandbox/multigrid/dec_test.py	2007-11-01 01:31:37 UTC (rev 3478)
@@ -1,10 +1,12 @@

from scipy import *
+from scipy.sparse import *
from pydec import *
-from pydec.multigrid import *
from pydec.multigrid.discrete_laplacian import boundary_hierarchy, discrete_laplacian_solver, hodge_solver

-from scipy.sandbox.multigrid import smoothed_aggregation_solver
+from scipy.sandbox.multigrid import smoothed_aggregation_solver,multigridtools,multilevel_solver
+from scipy.sandbox.multigrid.sa import sa_smoothed_prolongator
from scipy.sandbox.multigrid.utils import expand_into_blocks

@@ -13,8 +15,9 @@
-for i in range(3):
+for i in range(5):
mesh['vertices'],mesh['elements'] = loop_subdivision(mesh['vertices'],mesh['elements'])
cmplx = simplicial_complex(mesh['vertices'],mesh['elements'])

@@ -24,8 +27,34 @@
#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),dims=(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)),dims=(num_edges,num_basis*num_aggs)).tocsr()
+
+
+
+
+
def whitney_innerproduct_cache(cmplx,k):
h = hash(cmplx.vertices.tostring()) ^ hash(cmplx.simplices.tostring()) ^ hash(k)

@@ -73,47 +102,83 @@
Mi = whitney_innerproduct_cache(cmplx,i+1)
else:
Mi = regular_cube_innerproduct(cmplx,i+1)
+
+
+    dimension = mesh['vertices'].shape[1]

-    ##print "constructing solver"
-    ##ss = discrete_laplacian_solver(cochain_complex,len(cochain_complex)-i-1,innerproduct=Mi)
-    ##print ss
-    ##
-    ##print "solving"
-    ##x,res = ss.solve(b=zeros(ss.A.shape[0]),x0=rand(ss.A.shape[0]),return_residuals=True)
-
-    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
-    #A = B.T.tocsr() * Mi * B
-
-    constant_prolongators = [lvl[N - i].I for lvl in bh[:-1]]
-
-    if i == 0:
+    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:
-        #candidates = [ones(A.shape[0])]

-        #TODO test
-        candidates = []
-        for coord in range(mesh['vertices'].shape[1]):
-            candidates.append( bh[0][N-i+1].B * mesh['vertices'][:,coord] )
+        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'

-        K = len(candidates)
+        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:] ]

-        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':
+                ml = asa.solver
+            else:
+                raise ValuerError,'unknown method'
+
+        #ml = smoothed_aggregation_solver(A,candidates)

-
-    ml = smoothed_aggregation_solver(A,candidates,aggregation=constant_prolongators)
-    #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])
@@ -125,9 +190,11 @@
residuals.append(linalg.norm(b - A*x))
A.psolve = ml.psolve

+        from pydec import cg

+
residuals = array(residuals)/residuals[0]
avg_convergence_ratio = residuals[-1]**(1.0/len(residuals))
print "average convergence ratio",avg_convergence_ratio

Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-31 21:46:29 UTC (rev 3477)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2007-11-01 01:31:37 UTC (rev 3478)
@@ -24,14 +24,14 @@
O =  -ones(N)
return scipy.sparse.spdiags([D,O,O],[0,-1,1],N,N).tocoo().tocsr() #eliminate explicit zeros

-def poisson_problem2D(N):
+def poisson_problem2D(N,epsilon=1.0):
"""
Return a sparse CSR matrix for the 2d poisson problem
with standard 5-point finite difference stencil on a
square N-by-N grid.
"""
-    D = 4*ones(N*N)
-    T =  -ones(N*N)
+    D = (2 + 2*epsilon)*ones(N*N)
+    T =  -epsilon * ones(N*N)
O =  -ones(N*N)
T[N-1::N] = 0
return scipy.sparse.spdiags([D,O,T,T,O],[0,-N,-1,1,N],N*N,N*N).tocoo().tocsr() #eliminate explicit zeros
@@ -208,7 +208,9 @@
#use direct solver on coarsest level
#TODO reuse factors for efficiency?
coarse_x[:] = spsolve(self.As[-1],coarse_b).reshape(coarse_x.shape)
-            #coarse_x[:] = scipy.linalg.cg(self.As[-1],coarse_b,tol=1e-12)[0]
+            #coarse_x[:] = scipy.linalg.cg(self.As[-1],coarse_b,tol=1e-12)[0].reshape(coarse_x.shape)
+            #A_inv = asarray(scipy.linalg.pinv2(self.As[-1].todense()))
+            #coarse_x[:] = scipy.dot(A_inv,coarse_b)
#print "coarse residual norm",scipy.linalg.norm(coarse_b - self.As[-1]*coarse_x)
else:
self.__solve(lvl+1,coarse_x,coarse_b)
@@ -230,18 +232,17 @@
from scipy import *
candidates = None
blocks = None
-    #A = poisson_problem2D(100)
+    A = poisson_problem2D(40,1e-2)

-    #candidates = [ array(candidates[:,x]) for x in range(candidates.shape[1]) ]
-    blocks = arange(A.shape[0]/2).repeat(2)
+    #blocks = arange(A.shape[0]/2).repeat(2)

-    ml = smoothed_aggregation_solver(A,candidates,blocks=blocks,epsilon=0,max_coarse=100,max_levels=2)
+    ml = smoothed_aggregation_solver(A,candidates,blocks=blocks,epsilon=0.08,max_coarse=100,max_levels=10)
#ml = ruge_stuben_solver(A)

x = rand(A.shape[0])

Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py	2007-10-31 21:46:29 UTC (rev 3477)
+++ trunk/scipy/sandbox/multigrid/sa.py	2007-11-01 01:31:37 UTC (rev 3478)
@@ -182,7 +182,7 @@
P = T - (D_inv_A*T)

#S = (spidentity(A.shape[0]).tocsr() - D_inv_A) #TODO drop this?
-    #P = S * ( S * T)
+    #P = S *(S * ( S * T))

return P

@@ -201,8 +201,6 @@
T,coarse_candidates = sa_fit_candidates(AggOp,candidates)
#T = AggOp #TODO test

-    A_filtered = sa_filtered_matrix(A,epsilon,blocks) #use filtered matrix for anisotropic problems
-
P = sa_smoothed_prolongator(A,T,epsilon,omega,blocks)

if blocks is not None:

```