[Scipy-svn] r3854 - in trunk/scipy: sandbox/multigrid sparse

scipy-svn@scip... scipy-svn@scip...
Tue Jan 22 09:06:30 CST 2008


Author: wnbell
Date: 2008-01-22 09:06:22 -0600 (Tue, 22 Jan 2008)
New Revision: 3854

Modified:
   trunk/scipy/sandbox/multigrid/multilevel.py
   trunk/scipy/sandbox/multigrid/sa.py
   trunk/scipy/sparse/bsr.py
Log:
minor change to SA


Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2008-01-21 20:40:25 UTC (rev 3853)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2008-01-22 15:06:22 UTC (rev 3854)
@@ -102,35 +102,45 @@
 
     As = [A]
     Ps = []
+    Rs = []
 
     if aggregation is None:
         while len(As) < max_levels and A.shape[0] > max_coarse:
             P,B = sa_interpolation(A,B,epsilon*0.5**(len(As)-1),omega=omega)
+            R = P.T.asformat(P.format)
 
-            A = (P.T.asformat(P.format) * A) * P     #galerkin operator
+            A = R * A * P     #galerkin operator
 
             As.append(A)
+            Rs.append(R)
             Ps.append(P)
     else:
         #use user-defined aggregation
         for AggOp in aggregation:
             P,B = sa_interpolation(A,B,omega=omega,AggOp=AggOp)
+            R = P.T.asformat(P.format)
 
-            A = (P.T.tocsr() * A) * P     #galerkin operator
+            A = R * A * P     #galerkin operator
 
             As.append(A)
+            Rs.append(R)
             Ps.append(P)
 
-    return multilevel_solver(As,Ps,preprocess=pre,postprocess=post)
+    return multilevel_solver(As,Ps,Rs=Rs,preprocess=pre,postprocess=post)
 
 
 class multilevel_solver:
-    def __init__(self,As,Ps,preprocess=None,postprocess=None):
+    def __init__(self,As,Ps,Rs=None,preprocess=None,postprocess=None):
         self.As = As
         self.Ps = Ps
         self.preprocess = preprocess
         self.postprocess = postprocess
 
+        if Rs is None:
+            self.Rs = [P.T for P in self.Ps]
+        else:
+            self.Rs = Rs
+
     def __repr__(self):
         output = 'multilevel_solver\n'
         output += 'Number of Levels:     %d\n' % len(self.As)
@@ -202,7 +212,7 @@
 
         residual = b - A*x
 
-        coarse_b = self.Ps[lvl].T * residual
+        coarse_b = self.Rs[lvl] * residual
         coarse_x = zeros_like(coarse_b)
 
         if lvl == len(self.As) - 2:
@@ -212,7 +222,7 @@
             #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)
+            print "coarse residual norm",scipy.linalg.norm(coarse_b - self.As[-1]*coarse_x)
         else:
             self.__solve(lvl+1,coarse_x,coarse_b)
 

Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py	2008-01-21 20:40:25 UTC (rev 3853)
+++ trunk/scipy/sandbox/multigrid/sa.py	2008-01-22 15:06:22 UTC (rev 3854)
@@ -96,38 +96,38 @@
     else:
         sa_constant_interpolation(csr_matrix(A),epsilon)
 
-def sa_fit_candidates(AggOp,candidates,tol=1e-10):
+def sa_fit_candidates(AggOp,B,tol=1e-10):
     if not isspmatrix_csr(AggOp):
         raise TypeError,'expected csr_matrix for argument AggOp'
 
-    if candidates.dtype != 'float32':
-        candidates = asarray(candidates,dtype='float64')
+    if B.dtype != 'float32':
+        B = asarray(B,dtype='float64')
 
-    if len(candidates.shape) != 2:
+    if len(B.shape) != 2:
         raise ValueError,'expected rank 2 array for argument B'
 
-    if candidates.shape[0] % AggOp.shape[0] != 0:
+    if B.shape[0] % AggOp.shape[0] != 0:
         raise ValueError,'dimensions of AggOp %s and B %s are incompatible' % (AggOp.shape, B.shape)
     
 
-    K = candidates.shape[1] # number of near-nullspace candidates
-    blocksize = candidates.shape[0] / AggOp.shape[0]
+    K = B.shape[1] # number of near-nullspace candidates
+    blocksize = B.shape[0] / AggOp.shape[0]
 
     N_fine,N_coarse = AggOp.shape
 
-    R = zeros((N_coarse,K,K),dtype=candidates.dtype) #storage for coarse candidates
+    R = zeros((N_coarse,K,K), dtype=B.dtype) #storage for coarse candidates
 
     candidate_matrices = []
 
-    threshold = tol * abs(candidates).max()   # cutoff for small basis functions
-
     for i in range(K):
-        c = candidates[:,i]
+        c = B[:,i]
         c = c.reshape(-1,blocksize,1)[diff(AggOp.indptr) == 1]     # eliminate DOFs that aggregation misses
 
         X = bsr_matrix( (c, AggOp.indices, AggOp.indptr), \
                 shape=(blocksize*N_fine, N_coarse) )
 
+        col_thresholds = tol * bsr_matrix((X.data**2,X.indices,X.indptr),shape=X.shape).sum(axis=0).A.flatten() 
+
         #orthogonalize X against previous
         for j,A in enumerate(candidate_matrices):
             D_AtX = bsr_matrix((A.data*X.data,X.indices,X.indptr),shape=X.shape).sum(axis=0).A.flatten() #same as diagonal of A.T * X
@@ -135,9 +135,10 @@
             X.data -= scale_columns(A,D_AtX).data
 
         #normalize X
-        D_XtX = bsr_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   # set small basis functions to 0
+        col_norms = bsr_matrix((X.data**2,X.indices,X.indptr),shape=X.shape).sum(axis=0).A.flatten() #same as diagonal of X.T * X
+        mask = col_norms < col_thresholds   # set small basis functions to 0
+
+        col_norms = sqrt(col_norms)
         col_norms[mask] = 0
         R[:,i,i] = col_norms
         col_norms = 1.0/col_norms
@@ -186,7 +187,7 @@
 
     return P
 
-def sa_interpolation(A,candidates,epsilon=0.0,omega=4.0/3.0,AggOp=None):
+def sa_interpolation(A,B,epsilon=0.0,omega=4.0/3.0,AggOp=None):
 
     if not (isspmatrix_csr(A) or isspmatrix_bsr(A)):
         A = csr_matrix(A)
@@ -200,7 +201,7 @@
             raise ValueError,'incompatible aggregation operator'
 
 
-    T,coarse_candidates = sa_fit_candidates(AggOp,candidates)
+    T,coarse_candidates = sa_fit_candidates(AggOp,B)
     P = sa_smoothed_prolongator(A,T,epsilon,omega)
     return P,coarse_candidates
 

Modified: trunk/scipy/sparse/bsr.py
===================================================================
--- trunk/scipy/sparse/bsr.py	2008-01-21 20:40:25 UTC (rev 3853)
+++ trunk/scipy/sparse/bsr.py	2008-01-22 15:06:22 UTC (rev 3854)
@@ -166,7 +166,7 @@
             else:
                 self.shape = shape
 
-        self.check_format()
+        self.check_format(full_check=False)
 
     def check_format(self, full_check=True):
         """check whether the matrix format is valid



More information about the Scipy-svn mailing list