[Scipy-svn] r3431 - trunk/scipy/sandbox/multigrid
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Oct 10 13:12:39 EDT 2007
Author: wnbell
Date: 2007-10-10 12:12:37 -0500 (Wed, 10 Oct 2007)
New Revision: 3431
Added:
trunk/scipy/sandbox/multigrid/rs.py
trunk/scipy/sandbox/multigrid/sa.py
Modified:
trunk/scipy/sandbox/multigrid/adaptive.py
Log:
forgot to include rs.py and sa.py
Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-10 00:33:40 UTC (rev 3430)
+++ trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-10 17:12:37 UTC (rev 3431)
@@ -345,8 +345,8 @@
#A = io.mmread("/home/nathan/Desktop/BasisShift_W_EnergyMin_Luke/9pt-5x5.mtx").tocsr()
#A = A*A
-D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
-A = D * A * D
+#D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
+#A = D * A * D
#A = io.mmread("nos2.mtx").tocsr()
asa = adaptive_sa_solver(A,max_candidates=1)
#x = arange(A.shape[0]).astype('d') + 1
Added: trunk/scipy/sandbox/multigrid/rs.py
===================================================================
--- trunk/scipy/sandbox/multigrid/rs.py 2007-10-10 00:33:40 UTC (rev 3430)
+++ trunk/scipy/sandbox/multigrid/rs.py 2007-10-10 17:12:37 UTC (rev 3431)
@@ -0,0 +1,37 @@
+from scipy.sparse import csr_matrix,isspmatrix_csr
+
+import multigridtools
+
+__all__ = ['rs_strong_connections','rs_interpolation']
+
+def rs_strong_connections(A,theta):
+ """Return a strength of connection matrix using the method of Ruge and Stuben
+
+ An off-diagonal entry A[i.j] is a strong connection iff
+
+ -A[i,j] >= theta * max( -A[i,k] ) where k != i
+ """
+ if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
+
+ Sp,Sj,Sx = multigridtools.rs_strong_connections(A.shape[0],theta,A.indptr,A.indices,A.data)
+ return csr_matrix((Sx,Sj,Sp),dims=A.shape)
+
+
+def rs_interpolation(A,theta=0.25):
+ if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
+
+ S = rs_strong_connections(A,theta)
+
+ T = S.T.tocsr() #transpose S for efficient column access
+
+ Ip,Ij,Ix = multigridtools.rs_interpolation(A.shape[0],\
+ A.indptr,A.indices,A.data,\
+ S.indptr,S.indices,S.data,\
+ T.indptr,T.indices,T.data)
+
+ return csr_matrix((Ix,Ij,Ip))
+
+
+
+
+
Added: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py 2007-10-10 00:33:40 UTC (rev 3430)
+++ trunk/scipy/sandbox/multigrid/sa.py 2007-10-10 17:12:37 UTC (rev 3431)
@@ -0,0 +1,125 @@
+import scipy
+import numpy
+from numpy import array,arange,ones,zeros,sqrt,isinf,asarray,empty
+from scipy.sparse import csr_matrix,isspmatrix_csr
+
+from utils import diag_sparse,approximate_spectral_radius
+import multigridtools
+
+__all__ = ['sa_strong_connections','sa_constant_interpolation',
+ 'sa_interpolation','sa_fit_candidates']
+
+
+def sa_strong_connections(A,epsilon):
+ if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
+
+ Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
+
+ #D = diag_sparse(D)
+ #I,J,V = arange(A.shape[0]).repeat(diff(A.indptr)),A.indices,A.data #COO format for A
+ #diag_mask = (I == J)
+
+ return csr_matrix((Sx,Sj,Sp),A.shape)
+
+def sa_constant_interpolation(A,epsilon,blocks=None):
+ if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
+
+ #TODO handle epsilon = 0 case without creating strength of connection matrix?
+
+ if blocks is not None:
+ num_dofs = A.shape[0]
+ num_blocks = blocks.max()
+
+ if num_dofs != len(blocks):
+ raise ValueError,'improper block specification'
+
+ # for non-scalar problems, use pre-defined blocks in aggregation
+ # the strength of connection matrix is based on the Frobenius norms of the blocks
+
+ B = csr_matrix((ones(num_dofs),blocks,arange(num_dofs + 1)),dims=(num_dofs,num_blocks))
+ Block_Frob = B.T.tocsr() * csr_matrix((A.data**2,A.indices,A.indptr),dims=A.shape) * B #Frobenius norms of blocks entries of A
+
+ S = sa_strong_connections(Block_Frob,epsilon)
+
+ Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
+ Pj = Pj[blocks] #expand block aggregates into constituent dofs
+ Pp = B.indptr
+ Px = B.data
+ else:
+ S = sa_strong_connections(A,epsilon)
+
+ Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
+ Pp = numpy.arange(len(Pj)+1)
+ Px = numpy.ones(len(Pj))
+
+ return csr_matrix((Px,Pj,Pp))
+
+
+def sa_fit_candidates(AggOp,candidates):
+ K = len(candidates)
+
+ N_fine,N_coarse = AggOp.shape
+
+ 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))
+ N_fine = K*N_fine
+
+ #TODO convert this to list of coarse candidates
+ R = zeros((K*N_coarse,K)) #storage for coarse candidates
+
+ candidate_matrices = []
+ for i,c in enumerate(candidates):
+ #TODO permit incomplete AggOps here (for k-form problems) (other modifications necessary?)
+ X = csr_matrix((c.copy(),AggOp.indices,AggOp.indptr),dims=AggOp.shape)
+
+
+ #orthogonalize X against previous
+ for j,A in enumerate(candidate_matrices):
+ D_AtX = csr_matrix((A.data*X.data,X.indices,X.indptr),dims=X.shape).sum(axis=0).A.flatten() #same as diagonal of A.T * X
+ R[j::K,i] = D_AtX
+ X.data -= D_AtX[X.indices] * A.data
+
+
+ #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[i::K,i] = col_norms
+ col_norms = 1.0/col_norms
+ col_norms[isinf(col_norms)] = 0
+ X.data *= col_norms[X.indices]
+
+ candidate_matrices.append(X)
+
+ Q_indptr = K*AggOp.indptr
+ Q_indices = (K*AggOp.indices).repeat(K)
+ for i in range(K):
+ Q_indices[i::K] += i
+ Q_data = empty(N_fine * K)
+ for i,X in enumerate(candidate_matrices):
+ Q_data[i::K] = X.data
+ Q = csr_matrix((Q_data,Q_indices,Q_indptr),dims=(N_fine,K*N_coarse))
+
+ coarse_candidates = [array(R[:,i]) for i in range(K)]
+
+ return Q,coarse_candidates
+
+
+def sa_interpolation(A,candidates,epsilon,omega=4.0/3.0,blocks=None):
+ if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
+
+ AggOp = sa_constant_interpolation(A,epsilon=epsilon,blocks=blocks)
+ T,coarse_candidates = sa_fit_candidates(AggOp,candidates)
+
+ #TODO use filtered matrix here for anisotropic problems
+ A_filtered = A
+ 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)
+
+ P = T - (D_inv_A*T)
+
+ return P,coarse_candidates
+
+
+
More information about the Scipy-svn
mailing list