python/left_prm.py
changeset 0 0efde00f9229
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/python/left_prm.py	Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,92 @@
+#!/usr/bin/python
+
+import pylab as pl
+from rowjoin import rowjoin
+from polyfuncs import polsize
+from t1calc import t1calc
+from seshft import seshft
+from colsplit import colsplit
+from clcoef import clcoef
+
+def left_prm(N, degN, D, degD, job=1, gap=None):
+    """ """
+    if gap == None:
+        gap = 1.0e8
+    
+    F, degF = rowjoin(D, degD, N, degN)
+    Frows, Fbcols = polsize(F, degF)
+    Fcols = Fbcols * (degF+1)
+    T1 = pl.array([])
+    pr = pl.array([]) 
+    degT1, T1rows, shft = 0, 0, 0
+    S=F
+    sel = pl.ones((Frows,1))
+    T1bcols =1
+    abar = pl.arange(Fbcols, Frows)
+    if T1rows:
+        T1 = pl.hstack((T1, pl.zeros((T1rows, Frows))))
+    if T1.size == 0:
+        T1rows = 0
+    while T1.size == 0 or T1rows < (Frows - Fbcols):
+        if T1.size == 0:
+            T1rows = 0
+        Srows = Frows*T1bcols     
+        T1, T1rows, sel, pr = t1calc(S, Srows, T1, T1rows, sel, pr, Frows, Fbcols, abar, gap)
+        T1rows, T1cols = pl.atleast_2d(T1).shape
+        if T1.size == 0:
+            T1rows = 0
+        if T1rows < (Frows - Fbcols):
+            if T1.size != 0:
+                T1 = pl.hstack((T1, pl.zeros((T1rows,Frows)).squeeze()))
+            T1bcols += 1
+            degT1 += 1
+            shft += Fbcols
+            S = seshft(S,F,shft);
+            sel = pl.hstack((sel,sel[Srows-Frows:Srows+1]))
+            rowvec = pl.arange((T1bcols-1)*Frows+(Fbcols), (T1bcols * Frows))
+            abar = pl.hstack((abar, rowvec))
+
+    B, degB, A, degA = colsplit(T1, degT1, Fbcols, Frows-Fbcols)
+    B, degB = clcoef(B, degB)
+    B = -B;
+    A, degA = clcoef(A, degA)
+    if job == 1:
+        return B, degB, A, degA
+
+    if job == 2:
+        S = S[sel!=0,:]
+        redSrows, Scols = S.shape
+        C = pl.hstack(pl.eye(Fbcols), pl.zeros(Fbcols,Scols-Fbcols))
+        T2 = pl.dot(C, pl.inv(S))
+        T2 = makezero(T2,gap);
+        T2 = move(T2, pl.find(sel), Srows);
+        X, degX, Y, degY = colsplit(T2, degT1, Fbcols, Frows - Fbcols);
+        X, degX = clcoef(X, degX);
+        Y, degY = clcoef(Y, degY);
+    elif job == 3:
+        Y = S
+        degY = sel
+        X = degT1
+        degX = Fbcols
+
+    else:
+        print 'Message from left_prm:no legal job number specified'
+        exit()
+
+    return B, degB, A, degA, Y, degY, X, degX
+
+
+if __name__== "__main__":
+    D = pl.array([[1, 0, 0, 0, 0, 0],
+                  [0, 1, 0, 1, 0, 0],
+                  [0, 0, 1, 1, 1, 0]])
+
+    N = pl.array([[1, 0, 0],
+                  [0, 1, 0],
+                  [0, 0, 1]])
+
+    dD, dN = 1, 0
+
+    result = left_prm(N, dN, D, dD)
+    print result
+    pass