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

scipy-svn@scip... scipy-svn@scip...
Mon Oct 22 23:08:26 CDT 2007


Author: wnbell
Date: 2007-10-22 23:08:10 -0500 (Mon, 22 Oct 2007)
New Revision: 3458

Modified:
   trunk/scipy/sandbox/multigrid/adaptive.py
   trunk/scipy/sandbox/multigrid/multilevel.py
   trunk/scipy/sandbox/multigrid/relaxation.py
   trunk/scipy/sandbox/multigrid/utils.py
Log:
continued work on adaptiveSA


Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-22 20:13:01 UTC (rev 3457)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-23 04:08:10 UTC (rev 3458)
@@ -6,12 +6,66 @@
 from relaxation import gauss_seidel
 from multilevel import multilevel_solver
 from sa import sa_constant_interpolation,sa_fit_candidates
-from utils import approximate_spectral_radius,hstack_csr,vstack_csr
+from utils import approximate_spectral_radius,hstack_csr,vstack_csr,expand_into_blocks
 
 
+
+def augment_candidates(AggOp, old_Q, old_R, new_candidate, level):
+    K = old_R.shape[1]
     
+    #determine blocksizes
+    if level == 0:
+        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.shape[0]/new_bs[0],) + new_bs)
+    Q_data[:,:new_bs[0],:new_bs[1]] = old_Q.data.reshape((-1,) + old_bs)   #TODO BSR change 
 
+    # coarse candidates
+    R = zeros(((K+1)*AggOp.shape[1],K+1)) 
+    for i in range(K):
+        R[i::K+1,:K] = old_R[i::K,:] 
+
+    c = new_candidate.reshape(-1)[diff(AggOp.indptr) == 1]  #eliminate DOFs that aggregation misses
+    X = csr_matrix((c,AggOp.indices,AggOp.indptr),dims=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),dims=X.shape).sum(axis=0).A.flatten() #same as diagonal of A.T * X            
+        R[i::K+1,-1] = D_AtX
+        X.data -= D_AtX[X.indices] * old_c
+     
+    #normalize X
+    D_XtX = csr_matrix((X.data**2,X.indices,X.indptr),dims=X.shape).sum(axis=0).A.flatten() #same as diagonal of X.T * X            
+    col_norms = sqrt(D_XtX)
+    R[K::K+1,-1] = col_norms
+
+    col_norms = 1.0/col_norms
+    col_norms[isinf(col_norms)] = 0
+    X.data *= col_norms[X.indices]
+    last_column = Q_data[:,:,-1].reshape(-1)
+    last_column = X.data.reshape(-1)
+    Q_data = Q_data.reshape(-1)  #TODO BSR change
+
+    Q = csr_matrix((Q_data,Q_indices,Q_indptr),dims=(AggOp.shape[0],(K+1)*AggOp.shape[1]))
+
+    return Q,R
+    
+
+
+
 def orthonormalize_prolongator(P_l,x_l,W_l,W_m):
     """
      
@@ -72,26 +126,28 @@
     Construct multilevel hierarchy using Smoothed Aggregation
         Inputs:
           A  - matrix
-          Is - list of constant prolongators
+          Ps - list of constant prolongators
           B  - "candidate" basis function to be approximated
         Ouputs:
-          (As,Is,Ps) - tuple of lists
-                  - As - [A, Ps[0].T*A*Ps[0], Ps[1].T*A*Ps[1], ... ]
-                  - Is - smoothed prolongators 
-                  - Ps - tentative prolongators
+          (As,Ps,Ts) - tuple of lists
+                  - As - [A, Ts[0].T*A*Ts[0], Ts[1].T*A*Ts[1], ... ]
+                  - Ps - smoothed prolongators 
+                  - Ts - tentative prolongators
     """
     As = [A]
-    Is = []
     Ps = []
+    Ts = []
+    Bs = [B]
 
     for W in Ws:
         P,B = sa_fit_candidates(W,B)
         I   = smoothed_prolongator(P,A)  
         A   = I.T.tocsr() * A * I
         As.append(A)
-        Ps.append(P)
-        Is.append(I)
-    return As,Is,Ps
+        Ts.append(P)
+        Ps.append(I)
+        Bs.append(B)
+    return As,Ps,Ts,Bs
          
 def make_bridge(I,N):
     tail = I.indptr[-1].repeat(N - I.shape[0])
@@ -117,29 +173,30 @@
         
         Ws = AggOps
 
-        self.candidates = x
+        fine_candidates = x
 
         #create SA using x here
-        As,Is,Ps = sa_hierarchy(A,AggOps,self.candidates)
+        As,Ps,Ts,self.candidates = sa_hierarchy(A,AggOps,fine_candidates)
+        #self.candidates = [x]
 
         for i in range(max_candidates - 1):
-            #x = self.__develop_candidate(A,As,Is,Ps,Ws,AggOps,mu=mu)    
-            x = self.__develop_candidate(As,Is,Ps,AggOps,self.candidates,mu=mu)    
+            #x = self.__develop_candidate_OLD(A,As,Ps,Ts,Ws,AggOps,mu=mu)    
+            #As,Ps,Ts,Ws = self.__augment_cycle(A,As,Ts,Ws,AggOps,x)  
+            #self.candidates.append(x)
             
-            self.candidates = hstack((self.candidates,x))
-            
             #TODO which is faster?
-            #As,Is,Ps,Ws = self.__augment_cycle(A,As,Ps,Ws,AggOps,x)  
-            As,Is,Ps = sa_hierarchy(A,AggOps,self.candidates)
+            x = self.__develop_candidate(As,Ps,Ts,AggOps,self.candidates,mu=mu)    
+            fine_candidates = hstack((fine_candidates,x))
+            As,Ps,Ts,self.candidates = sa_hierarchy(A,AggOps,fine_candidates)
 
-        self.Ps = Ps 
-        self.solver = multilevel_solver(As,Is)
+        self.Ts = Ts 
+        self.solver = multilevel_solver(As,Ps)
         self.AggOps = AggOps
                 
 
     def __initialization_stage(self,A,blocks,max_levels,max_coarse,mu,epsilon):
         AggOps = []
-        Is     = []
+        Ps     = []
 
         # aSA parameters 
         # mu      - number of test relaxation iterations
@@ -166,7 +223,7 @@
            
             print 'A_l.shape',A_l.shape
             AggOps.append(W_l)
-            Is.append(I_l)
+            Ps.append(I_l)
             As.append(A_l)
 
             if A_l.shape <= max_coarse:  break
@@ -183,7 +240,7 @@
                         x = x_hat  #need to restore x
         
         #update fine-level candidate
-        for A_l,I in reversed(zip(As[1:],Is)):
+        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
@@ -191,93 +248,113 @@
         return x,AggOps  #first candidate,aggregation
 
 
-    def __develop_candidate(self,As,Is,Ps,AggOps,candidates,mu):
-        #scipy.random.seed(0)  #TEST
-        x = scipy.rand(A.shape[0])
+    def __develop_candidate(self,As,Ps,Ts,AggOps,candidates,mu):
+        A = As[0]
+        
+        x = A*scipy.rand(A.shape[0],1)
         b = zeros_like(x)
 
-        solver = multilevel_solver(As,Is)
-
-        x = solver.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
 
-        #TODO augment candiates each time, then call fit_candidates?
+        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),dims=(K*P.shape[0]/(K-1),P.shape[1]))
 
-        A_l,P_l,W_l,x_l = As[0],Ps[0],Ws[0],x
-        
-        temp_Is = []
-        for i in range(len(As) - 2):
-            P_l_new, x_m, W_l_new, W_m_new = orthonormalize_prolongator(P_l, x_l, W_l, AggOps[i+1])    
+        for i in range(len(As) - 2):            
+            #TODO test augment_candidates against fit candidates
+            if i == 0:
+                temp = hstack((candidates[0],x))
+            else:
+                K = candidates[i].shape[1] 
+                temp = zeros((x.shape[0]/(K+1),K+1,K + 1))
+                temp[:,:-1,:-1] = candidates[i].reshape(-1,K,K)
+                temp[:,:,-1] = x.reshape(-1,K+1,1)
+                temp = temp.reshape(-1,K+1)
+            T_,R_ = sa_fit_candidates(AggOps[i],temp)
+            #print "T - T_",(T - T_).data.max()
+            #assert((T - T_).data.max() < 1e-10)
+            #assert((R - R_).data.max() < 1e-10)
+            T,R = T_,R_
+            #TODO end test
+
+            #T,R = augment_candidates(AggOps[i], Ts[i], candidates[i+1], x, i)
+            P = 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]) 
  
-            I_l_new = smoothed_prolongator(P_l_new,A_l)
-            A_m_new = I_l_new.T.tocsr() * A_l * I_l_new
-            bridge = make_bridge(Is[i+1],A_m_new.shape[0]) 
+            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)
  
-            temp_solver = multilevel_solver( [A_m_new] + As[i+2:], [bridge] + Is[i+2:] )
- 
-            for n in range(mu):
-                x_m = temp_solver.solve(zeros_like(x_m), x0=x_m, tol=1e-8, maxiter=1)
- 
-            temp_Is.append(I_l_new)
- 
-            W_l = vstack_csr(Ws[i+1],W_m_new)  #prepare for next iteration
-            A_l = A_m_new
-            x_l = x_m
-            P_l = make_bridge(Ps[i+1],A_m_new.shape[0])
- 
-        x = x_l
-        for I in reversed(temp_Is):
-            x = I*x
+        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 P in reversed(temp_Ps):
+        #    x = P*x
+        
         return x
 
 
-##    def __develop_candidate(self,A,As,Is,Ps,Ws,AggOps,mu):
+##    def __develop_candidate_OLD(self,A,As,Ps,Ts,Ws,AggOps,mu):
 ##        #scipy.random.seed(0)  #TEST
 ##        x = scipy.rand(A.shape[0])
 ##        b = zeros_like(x)
 ##
-##        solver = multilevel_solver(As,Is)
+##        solver = multilevel_solver(As,Ps)
 ##
 ##        x = solver.solve(b, x0=x, tol=1e-10, maxiter=mu)
 ##    
 ##        #TEST FOR CONVERGENCE HERE
 ## 
-##        A_l,P_l,W_l,x_l = As[0],Ps[0],Ws[0],x
+##        A_l,P_l,W_l,x_l = As[0],Ts[0],Ws[0],x
 ##        
-##        temp_Is = []
+##        temp_Ps = []
 ##        for i in range(len(As) - 2):
 ##            P_l_new, x_m, W_l_new, W_m_new = orthonormalize_prolongator(P_l, x_l, W_l, AggOps[i+1])    
 ## 
 ##            I_l_new = smoothed_prolongator(P_l_new,A_l)
 ##            A_m_new = I_l_new.T.tocsr() * A_l * I_l_new
-##            bridge = make_bridge(Is[i+1],A_m_new.shape[0]) 
+##            bridge = make_bridge(Ps[i+1],A_m_new.shape[0]) 
 ## 
-##            temp_solver = multilevel_solver( [A_m_new] + As[i+2:], [bridge] + Is[i+2:] )
+##            temp_solver = multilevel_solver( [A_m_new] + As[i+2:], [bridge] + Ps[i+2:] )
 ## 
 ##            for n in range(mu):
 ##                x_m = temp_solver.solve(zeros_like(x_m), x0=x_m, tol=1e-8, maxiter=1)
 ## 
-##            temp_Is.append(I_l_new)
+##            temp_Ps.append(I_l_new)
 ## 
 ##            W_l = vstack_csr(Ws[i+1],W_m_new)  #prepare for next iteration
 ##            A_l = A_m_new
 ##            x_l = x_m
-##            P_l = make_bridge(Ps[i+1],A_m_new.shape[0])
+##            P_l = make_bridge(Ts[i+1],A_m_new.shape[0])
 ## 
 ##        x = x_l
-##        for I in reversed(temp_Is):
+##        for I in reversed(temp_Ps):
 ##            x = I*x
 ##        
 ##        return x
            
 
-    def __augment_cycle(self,A,As,Ps,Ws,AggOps,x):
+    def __augment_cycle(self,A,As,Ts,Ws,AggOps,x):
         #make a new cycle using the new candidate
-        A_l,P_l,W_l,x_l = As[0],Ps[0],AggOps[0],x            
+        A_l,P_l,W_l,x_l = As[0],Ts[0],AggOps[0],x            
         
-        new_As,new_Is,new_Ps,new_Ws = [A],[],[],[AggOps[0]]
+        new_As,new_Ps,new_Ts,new_Ws = [A],[],[],[AggOps[0]]
 
         for i in range(len(As) - 2):
             P_l_new, x_m, W_l_new, W_m_new = orthonormalize_prolongator(P_l, x_l, W_l, AggOps[i+1])
@@ -288,24 +365,24 @@
 
             new_As.append(A_m_new)
             new_Ws.append(W_m_new)
-            new_Is.append(I_l_new)
-            new_Ps.append(P_l_new)
+            new_Ps.append(I_l_new)
+            new_Ts.append(P_l_new)
 
             #prepare for next iteration
             W_l = W_m_new 
             A_l = A_m_new
             x_l = x_m
-            P_l = make_bridge(Ps[i+1],A_m_new.shape[0]) 
+            P_l = make_bridge(Ts[i+1],A_m_new.shape[0]) 
         
         P_l_new, x_m, W_l_new, W_m_new = orthonormalize_prolongator(P_l, x_l, W_l, csr_matrix((P_l.shape[1],1)))
         I_l_new = smoothed_prolongator(P_l_new,A_l)
         A_m_new = I_l_new.T.tocsr() * A_l * I_l_new
 
         new_As.append(A_m_new)
-        new_Is.append(I_l_new)
-        new_Ps.append(P_l_new)
+        new_Ps.append(I_l_new)
+        new_Ts.append(P_l_new)
 
-        return new_As,new_Is,new_Ps,new_Ws
+        return new_As,new_Ps,new_Ts,new_Ws
 
 
 
@@ -325,10 +402,10 @@
 #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)
+A = io.mmread("tests/sample_data/elas30_A.mtx").tocsr()
+blocks = arange(A.shape[0]/2).repeat(2)
 
-asa = adaptive_sa_solver(A,max_candidates=1,mu=10)
+asa = adaptive_sa_solver(A,max_candidates=4,mu=5)
 scipy.random.seed(0)  #make tests repeatable
 x = rand(A.shape[0])
 #b = A*rand(A.shape[0])
@@ -356,7 +433,7 @@
 
 print "constant Rayleigh quotient",dot(ones(A.shape[0]),A*ones(A.shape[0]))/float(A.shape[0])
 
-for c in asa.candidates.T:
+for c in asa.candidates[0].T:
     print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
 
 
@@ -364,6 +441,19 @@
 ##W = asa.AggOps[0]*asa.AggOps[1]
 ##pcolor((W * rand(W.shape[1])).reshape((200,200)))
 
+def plot2d_arrows(x):
+    from pylab import quiver
+    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]
+
+    quiver(X,Y,dX,dY)
+
 def plot2d(x):
     from pylab import pcolor
     pcolor(x.reshape(sqrt(len(x)),sqrt(len(x))))

Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-22 20:13:01 UTC (rev 3457)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-23 04:08:10 UTC (rev 3458)
@@ -168,7 +168,7 @@
         if x0 is None:
             x = zeros_like(b)
         else:
-            x = array(x0)
+            x = array(x0) #copy
 
 
         #TODO change use of tol (relative tolerance) to agree with other iterative solvers
@@ -193,19 +193,21 @@
         A = self.As[lvl]
         
         if len(self.As) == 1:
-            x[:] = spsolve(A,b)
+            #TODO make spsolve preserve dimensions
+            x[:] = spsolve(A,b).reshape(x.shape)
             return 
 
         self.presmoother(A,x,b)
 
         residual = b - A*x                                
 
-        coarse_x = zeros((self.As[lvl+1].shape[0]))
         coarse_b = self.Ps[lvl].T * residual
+        coarse_x = zeros_like(coarse_b)
         
         if lvl == len(self.As) - 2:
             #use direct solver on coarsest level
-            coarse_x[:] = spsolve(self.As[-1],coarse_b)  #TODO reuse factors for efficiency?
+            #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]
             #print "coarse residual norm",scipy.linalg.norm(coarse_b - self.As[-1]*coarse_x)
         else:   

Modified: trunk/scipy/sandbox/multigrid/relaxation.py
===================================================================
--- trunk/scipy/sandbox/multigrid/relaxation.py	2007-10-22 20:13:01 UTC (rev 3457)
+++ trunk/scipy/sandbox/multigrid/relaxation.py	2007-10-23 04:08:10 UTC (rev 3458)
@@ -77,6 +77,9 @@
          iterations - number of iterations to perform (default: 1)
          omega      - damping parameter (default: 1.0)
     """ 
+    x = asarray(x).reshape(-1)
+    b = asarray(b).reshape(-1)
+
     sweep = slice(None)
     (row_start,row_stop,row_step) = sweep.indices(A.shape[0])
     

Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py	2007-10-22 20:13:01 UTC (rev 3457)
+++ trunk/scipy/sandbox/multigrid/utils.py	2007-10-23 04:08:10 UTC (rev 3458)
@@ -65,8 +65,8 @@
     if not isspmatrix(A) or not isspmatrix(B):
         raise TypeError,'expected sparse matrix'
     
-    if A.shape[0] != B.shape[0]:
-        raise ValueError,'row dimensions must agree'
+    if A.shape[1] != B.shape[1]:
+        raise ValueError,'column dimensions must agree'
 
     A = A.tocoo()
     B = B.tocoo()
@@ -94,9 +94,6 @@
     #TODO EXPLAIN MORE
     #TODO use spkron instead, time for compairson
 
-    if n is None:
-        n = m
-
     if m == 1 and n == 1:
         return A #nothing to do
 



More information about the Scipy-svn mailing list