[Scipy-svn] r3445 - trunk/scipy/sandbox/multigrid

scipy-svn@scip... scipy-svn@scip...
Thu Oct 18 03:06:14 CDT 2007


Author: wnbell
Date: 2007-10-18 03:06:05 -0500 (Thu, 18 Oct 2007)
New Revision: 3445

Added:
   trunk/scipy/sandbox/multigrid/dec_test.py
Modified:
   trunk/scipy/sandbox/multigrid/adaptive.py
   trunk/scipy/sandbox/multigrid/multilevel.py
   trunk/scipy/sandbox/multigrid/sa.py
   trunk/scipy/sandbox/multigrid/utils.py
Log:
fixed bug in expansion of AggOp
added dec example, incomplete 



Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-17 20:29:14 UTC (rev 3444)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-18 08:06:05 UTC (rev 3445)
@@ -294,7 +294,7 @@
 
 print "solving"
 if True:
-    x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=10,tol=1e-12,return_residuals=True)
+    x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=25,tol=1e-8,return_residuals=True)
 else:
     residuals = []
     def add_resid(x):

Added: trunk/scipy/sandbox/multigrid/dec_test.py
===================================================================
--- trunk/scipy/sandbox/multigrid/dec_test.py	2007-10-17 20:29:14 UTC (rev 3444)
+++ trunk/scipy/sandbox/multigrid/dec_test.py	2007-10-18 08:06:05 UTC (rev 3445)
@@ -0,0 +1,174 @@
+
+from scipy 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.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')
+for i in range(3):
+    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 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)
+
+    ##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:
+        candidates = None
+    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] )
+
+        K = len(candidates)
+
+        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)
+    #ml = smoothed_aggregation_solver(A,candidates)
+
+    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
+        x_sol = linalg.cg(A,b,x0=x,maxiter=30,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
+
+
+   
+
+##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
+##

Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-17 20:29:14 UTC (rev 3444)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-18 08:06:05 UTC (rev 3445)
@@ -87,7 +87,7 @@
     else:
         #use user-defined aggregation
         for AggOp in aggregation:
-            P,candidates = sa_interpolation(A,candidates,omega=omega,AggOp=AggOp)
+            P,candidates,blocks = sa_interpolation(A,candidates,omega=omega,AggOp=AggOp)
 
             A = (P.T.tocsr() * A) * P     #galerkin operator
 
@@ -210,11 +210,11 @@
     #ml = ruge_stuben_solver(A)
 
     x = rand(A.shape[0])
-    #b = zeros_like(x)
-    b = A*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)
+        x_sol,residuals = ml.solve(b,x0=x,maxiter=30,tol=1e-8,return_residuals=True)
     else:
         residuals = []
         def add_resid(x):

Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py	2007-10-17 20:29:14 UTC (rev 3444)
+++ trunk/scipy/sandbox/multigrid/sa.py	2007-10-18 08:06:05 UTC (rev 3445)
@@ -2,28 +2,53 @@
 import numpy
 from numpy import array,arange,ones,zeros,sqrt,isinf,asarray,empty,diff,\
                   ascontiguousarray
-from scipy.sparse import csr_matrix,isspmatrix_csr
+from scipy.sparse import csr_matrix,isspmatrix_csr,spidentity
 
-from utils import diag_sparse,approximate_spectral_radius
+from utils import diag_sparse,approximate_spectral_radius,expand_into_blocks
 import multigridtools
 
 __all__ = ['sa_filtered_matrix','sa_strong_connections','sa_constant_interpolation',
            'sa_interpolation','sa_smoothed_prolongator','sa_fit_candidates']
 
 
+##    nnz = A.nnz
+##
+##    indptr  = A.indptr
+##    indices = A.indices
+##    data    = A.data
+##
+##    if n != 1:
+##        # expand horizontally
+##        indptr = n*A.indptr
+##        indices = (n*indices).repeat(n) + tile(arange(n),nnz) 
+##
+##    if m != 1:
+##        #expand vertically
+##        indptr  = concatenate( (array([0]), cumsum(diff(A).repeat(m))) )
+##        indices = indices.repeat(m)
+##
+##    if m != 1 or n != 1:
+##        data = A.data.repeat(m*n)
+
+
 def sa_filtered_matrix(A,epsilon,blocks=None):
+    """The filtered matrix is obtained from A by lumping all weak off-diagonal 
+    entries onto the diagonal.  Weak off-diagonals are determined by
+    the standard strength of connection measure using the parameter epsilon.
+
+    In the case epsilon = 0.0, (i.e. no weak connections) A is returned.
+    """
+
     if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
 
     if epsilon == 0:
         A_filtered = A
-
     else:
         if blocks is None:
             Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
             A_filtered = csr_matrix((Sx,Sj,Sp),A.shape)
         else:
             A_filtered = A  #TODO subtract weak blocks from diagonal blocks?
-
 ##            num_dofs   = A.shape[0]
 ##            num_blocks = blocks.max() + 1
 ##            
@@ -50,7 +75,6 @@
 
     return A_filtered          
 
-
 def sa_strong_connections(A,epsilon):
     if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
     
@@ -101,7 +125,9 @@
 
     if K > 1 and len(candidates[0]) == K*N_fine:
         #see if fine space has been expanded (all levels except for first)
-        AggOp = csr_matrix((AggOp.data.repeat(K),AggOp.indices.repeat(K),arange(K*N_fine + 1)),dims=(K*N_fine,N_coarse))
+        #TODO fix this and add unittests
+        #AggOp = csr_matrix((AggOp.data.repeat(K),AggOp.indices.repeat(K),arange(K*N_fine + 1)),dims=(K*N_fine,N_coarse))
+        AggOp = expand_into_blocks(AggOp,K,1).tocsr()
         N_fine = K*N_fine
 
     #TODO convert this to list of coarse candidates
@@ -142,6 +168,20 @@
     return Q,coarse_candidates
     
 def sa_smoothed_prolongator(A,T,epsilon,omega,blocks=None):
+    """For a given matrix A and tentative prolongator T return the 
+    smoothed prolongator P
+
+        P = (I - omega/rho(S) S) * T
+
+    where S is a Jacobi smoothing operator defined as follows:
+
+        omega      - damping parameter
+        rho(S)     - spectral radius of S (estimated)
+        S          - inv(diag(A_filtered)) * A_filtered   (Jacobi smoother)
+        A_filtered - sa_filtered_matrix(A,epsilon)
+    """
+     
+
     A_filtered = sa_filtered_matrix(A,epsilon,blocks) #use filtered matrix for anisotropic problems
 
     D_inv    = diag_sparse(1.0/diag_sparse(A_filtered))       
@@ -150,6 +190,9 @@
 
     # smooth tentative prolongator T
     P = T - (D_inv_A*T)
+    
+    #S = (spidentity(A.shape[0]).tocsr() - D_inv_A) #TODO drop this?
+    #P = S * ( S * T)
 
     return P
 
@@ -163,20 +206,15 @@
             raise TypeError,'expected csr_matrix for argument AggOp'
         if A.shape[1] != AggOp.shape[0]:
             raise ValueError,'incompatible aggregation operator'
+        
 
     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)
 
-##    D_inv    = diag_sparse(1.0/diag_sparse(A_filtered))       
-##    D_inv_A  = D_inv * A_filtered
-##    D_inv_A *= omega/approximate_spectral_radius(D_inv_A)
-##
-##    # smooth tentative prolongator T
-##    P = T - (D_inv_A*T)
-
     if blocks is not None:
         blocks = arange(AggOp.shape[1]).repeat(len(candidates))
           

Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py	2007-10-17 20:29:14 UTC (rev 3444)
+++ trunk/scipy/sandbox/multigrid/utils.py	2007-10-18 08:06:05 UTC (rev 3445)
@@ -1,9 +1,9 @@
 __all__ =['approximate_spectral_radius','infinity_norm','diag_sparse',
-          'hstack_csr','vstack_csr']
+          'hstack_csr','vstack_csr','expand_into_blocks']
 
 import numpy
 import scipy
-from numpy import ravel,arange,concatenate
+from numpy import ravel,arange,concatenate,tile
 from scipy.linalg import norm
 from scipy.sparse import isspmatrix,isspmatrix_csr,isspmatrix_csc, \
                         csr_matrix,csc_matrix,extract_diagonal, \
@@ -67,3 +67,47 @@
     V = concatenate((A.data,B.data))
     return coo_matrix((V,(I,J)),dims=(A.shape[0]+B.shape[0],A.shape[1])).tocsr()
 
+
+def expand_into_blocks(A,m,n):
+    """Expand each element in a sparse matrix A into an m-by-n block.  
+                
+          Example: 
+          >>> A.todense()
+          matrix([[ 1.,  2.],
+                  [ 4.,  5.]])
+          
+          >>> expand_into_blocks(A,2,2).todense()
+          matrix([[ 1.,  1.,  2.,  2.],
+                  [ 1.,  1.,  2.,  2.],
+                  [ 4.,  4.,  5.,  5.],
+                  [ 4.,  4.,  5.,  5.]])
+              
+    """
+    #TODO EXPLAIN MORE
+
+    if n is None:
+        n = m
+
+    if m == 1 and n == 1:
+        return A #nothing to do
+
+    A = A.tocoo()
+
+    # expand 1x1 -> mxn
+    row  = ( m*A.row ).repeat(m*n).reshape(-1,m,n)
+    col  = ( n*A.col ).repeat(m*n).reshape(-1,m,n)
+
+    # increment indices
+    row += tile(arange(m).reshape(-1,1),(1,n))
+    col += tile(arange(n).reshape(1,-1),(m,1))
+
+    # flatten
+    row = row.reshape(-1)
+    col = col.reshape(-1)
+
+    data = A.data.repeat(m*n)
+
+    return coo_matrix((data,(row,col)),dims=(m*A.shape[0],n*A.shape[1]))
+
+
+



More information about the Scipy-svn mailing list