def dampnewton(x,F,DF,q=0.5,tol=1e-10):
    cvg = []
    lup = lu_factor(DF(x))
    s = lu_solve(lup,F(x))
    xn = x-s
    lam = 1
    st = lu_solve(lup,F(xn)) # simplified Newton
    while norm(st) > tol*norm(xn):
        while norm(st) > (1-lam*0.5)*norm(s):
            lam *= 0.5
            if lam < 1e-10:
                cvg = -1
                print 'Failure of convergence'
                return x, cvg
            xn = x-lam*s
            st = lu_solve(lup,F(xn)) # simplified Newton
        cvg += [[lam, norm(xn), norm(F(xn))]]
        x = xn
        lup = lu_factor(DF(x))
        s = lu_solve(lup,F(x))
        lam = min(lam/q, 1.) # Wozu dieser Test?
        xn = x-lam*s
        st = lu_solve(lup,F(xn)) # simplified Newton
    x = xn
    return x, array(cvg)

    def __init__(self, N, a0, alfa, beta, quad="GL", solver="scipy"):
        self.quad = quad
        self.solver = solver
        k = arange(N)
        self.S = S = SBBmat(k)
        self.B = B = BBBmat(k, self.quad)
        self.A = A = ABBmat(k)
        self.a0 = a0
        self.alfa = alfa
        self.beta = beta
        if not solver == "scipy":
            sii, siu, siuu = S.dd, S.ud[0], S.ud[1]
            ail, aii, aiu = A.ld, A.dd, A.ud
            bill, bil, bii, biu, biuu = B.lld, B.ld, B.dd, B.ud, B.uud
            M = sii[::2].shape[0]
        
        if hasattr(beta, "__len__"):
            Ny, Nz = beta.shape
            if solver == "scipy":
                self.Le = Le = []
                self.Lo = Lo = []
                for i in range(Ny):
                    Lej = []
                    Loj = []
                    for j in range(Nz):
                        AA = a0*S.diags().toarray() + alfa[i, j]*A.diags().toarray() + beta[i, j]*B.diags().toarray()
                        Ae = AA[::2, ::2]
                        Ao = AA[1::2, 1::2]
                        Lej.append(lu_factor(Ae))
                        Loj.append(lu_factor(Ao))
                    Le.append(Lej)
                    Lo.append(Loj)
            else:
                self.u0 = zeros((2, M, Ny, Nz))
                self.u1 = zeros((2, M, Ny, Nz))
                self.u2 = zeros((2, M, Ny, Nz))
                self.l0 = zeros((2, M, Ny, Nz))
                self.l1 = zeros((2, M, Ny, Nz))
                self.ak = zeros((2, M, Ny, Nz))
                self.bk = zeros((2, M, Ny, Nz))
                SFTc.LU_Biharmonic_3D(a0, alfa, beta, sii, siu, siuu, ail, aii, aiu, bill, bil, bii, biu, biuu, self.u0, self.u1, self.u2, self.l0, self.l1)
                SFTc.Biharmonic_factor_pr_3D(self.ak, self.bk, self.l0, self.l1)

        else:
            if solver == "scipy":
                AA = a0*S.diags().toarray() + alfa*A.diags().toarray() + beta*B.diags().toarray()
                Ae = AA[::2, ::2]
                Ao = AA[1::2, 1::2]
                self.Le = lu_factor(Ae)
                self.Lo = lu_factor(Ao)
            else:
                self.u0 = zeros((2, M))
                self.u1 = zeros((2, M))
                self.u2 = zeros((2, M))
                self.l0 = zeros((2, M))
                self.l1 = zeros((2, M))
                self.ak = zeros((2, M))
                self.bk = zeros((2, M))
                SFTc.LU_Biharmonic_1D(a0, alfa, beta, sii, siu, siuu, ail, aii, aiu, bill, bil, bii, biu, biuu, self.u0, self.u1, self.u2, self.l0, self.l1)
                SFTc.Biharmonic_factor_pr(self.ak, self.bk, self.l0, self.l1)

    def solve_nonlinear(self, inputs, outputs):
        """
        Use numpy to solve Ax=b for x.

        Parameters
        ----------
        inputs : Vector
            unscaled, dimensional input variables read via inputs[key]
        outputs : Vector
            unscaled, dimensional output variables read via outputs[key]
        """
        vec_size = self.options['vec_size']
        vec_size_A = self.vec_size_A

        # lu factorization for use with solve_linear
        self._lup = []
        if vec_size > 1:
            for j in range(vec_size_A):
                lhs = inputs['A'][j] if vec_size_A > 1 else inputs['A']
                self._lup.append(linalg.lu_factor(lhs))

            for j in range(vec_size):
                idx = j if vec_size_A > 1 else 0
                outputs['x'][j] = linalg.lu_solve(self._lup[idx], inputs['b'][j])
        else:
            self._lup = linalg.lu_factor(inputs['A'])
            outputs['x'] = linalg.lu_solve(self._lup, inputs['b'])

def build_np_models(kernel, trans_samples, ter_samples, ter_rew_samples, lamb):
    Xa, Ra, Xpa = zip(*trans_samples)
    Xa_term, Ra_term = zip(*ter_samples)
    Xa = [ np.vstack((xa, xa_term)) if xa_term.size > 0 else xa for xa, xa_term in izip(Xa, Xa_term) ]
    Ra = [ np.hstack((ra, ra_term)) if ra_term.size > 0 else ra for ra, ra_term in izip(Ra, Ra_term) ] 
    
    k = len(trans_samples) 
    
    # build the K_a,b matrices
    Kab = dict()
    for a,b in product(xrange(k), xrange(k)):
        if Xa_term[b].size > 0:
            Kab[(a,b)] = np.hstack((kernel(Xa[a], Xpa[b]), 
                                    np.zeros((Xa[a].shape[0], Xa_term[b].shape[0]))))
        else:
            Kab[(a,b)] = kernel(Xa[a], Xpa[b])
        
    
    # build the K_a, D_a matrices
    Ka = [kernel(Xa[i], Xa[i])  for i in xrange(k)]
    Dainv = [Ka[i] + lamb*scipy.eye(*Ka[i].shape) for i in xrange(k)]
    Da = [lu_factor(Dainv[i], overwrite_a = False) for i in xrange(k)]
        
    # build K_ter matrix
    Kterma = [ np.hstack((kernel(ter_rew_samples[0], Xpa[i]),
                          np.zeros((ter_rew_samples[0].shape[0], Xa_term[i].shape[0])))) if Xa_term[i].size > 0
                else kernel(ter_rew_samples[0], Xpa[i]) for i in xrange(k)]
    K_ter = kernel(ter_rew_samples[0], ter_rew_samples[0])
    D_ter = lu_factor(K_ter + lamb*scipy.eye(*K_ter.shape), overwrite_a = True)
    R_ter = ter_rew_samples[1]
    
    return kernel, Kab, Da, Dainv, Ra, Kterma, D_ter, R_ter, Xa

    def __call__(self, A, shift=0):
        from scipy.linalg import lu_factor, lu_solve
        nrows = A.mshape[0]
        nblock = A.blockshape[0] #assume a square block!

        assert A.bandwidth==3, "Matrix bust be tridiagonal block matrix"

        #Overwrite current matrix?
        if self.overwrite:
            self.Alu = A
            
        #Create new internal matrix if none
        elif not hasattr(self, 'Alu'):
            logging.debug( "Creating new internal LU matrix" )
            self.Alu=BlockArray(A.mshape, A.blockshape, dtype=A.dtype)

        #Create new internal matrix if A is different shape
        elif (self.Alu.mshape<>A.mshape) or (self.Alu.blockshape<>A.blockshape):
            logging.debug( "Internal LU matrix incorrect shape; creating new one" )
            del self.Alu
            self.Alu = BlockArray(A.mshape, A.blockshape, dtype=A.dtype)

        #Vector for pivots
        self.shift = shift
        self.pivots = zeros((nrows, nblock), dtype=A.dtype)
        piv=0
                
        Ishift = shift*eye(nblock)
        self.Alu.set_raw_block(0,1,A.get_raw_block(0, 1) - Ishift)
        bnew = self.Alu.get_raw_block(0, 1)
        for row in range(0,nrows-1):
            a = A.get_raw_block(row+1, 0)
            b = bnew
            c = A.get_raw_block(row, 2)
            
            b, self.pivots[row] = lu_factor(b)
            # anew = a inv(b)
            anew = lu_solve((b,self.pivots[row]), a.T, trans=1).T
            bnew = A.get_raw_block(row+1, 1) - Ishift - dot(anew, c) 

            self.Alu.set_raw_block(row,1,b)                 #blu
            self.Alu.set_raw_block(row+1,0,anew)    #b anew = a

            if not self.overwrite:                      #Copy over block c, if not overwriting
                                        self.Alu.set_raw_block(row,2,c)

        #Final LU decomp of last block
        b, self.pivots[nrows-1] = lu_factor(bnew)
        self.Alu.set_raw_block(nrows-1,1,b)

    def linearize(self, params, unknowns, resids):
        """ Jacobian for Sellar discipline 1."""

        z1 = params['z'][0]
        z2 = params['z'][1]
        x1 = params['x']
        y1 = unknowns['y1']
        y2 = unknowns['y2']

        J = {}

        J['y1', 'y2'] = -0.2
        J['y1', 'z'] = np.array([[2.0*z1, 1.0]])
        J['y1', 'x'] = 1.0
        J['y2', 'y1'] = .5*y1**-.5
        J['y2', 'z'] = np.array([[1.0, 1.0]])
        J['y2', 'y2'] = -1.0

        dRdy = np.zeros((2, 2))
        dRdy[0, 1] = J['y1', 'y2']
        dRdy[0, 0] = 1.0
        dRdy[1, 0] = J['y2', 'y1']
        dRdy[1, 1] = J['y2', 'y2']

        # lu factorization for use with solve_linear
        self.lup = linalg.lu_factor(dRdy)

        return J

def factiz(K):
    """
    Helper function to behave the same way scipy.sparse.factorized does, but
    for dense matrices.
    """
    luf = lu_factor(K)
    return lambda x: matrix(lu_solve(luf, x))

def invpower(e,A,B=None):
	if B is None:
		B = eye(len(A))
	K = A - B*e
	G = lu_factor(K)
	x = ones([len(A),1],complex)/len(A)
	iter = 0
	error = 10
	#for i in range(8):
	while error > 1e-8 and iter < 20:
		try:
			x = lu_solve(G,x)
		    
		except:
			print 'LU Exception'
			x = solve(K,x)
		    
		x = x/norm(x,ord=inf)
		error = norm(dot(A,x)-e*dot(B,x))
		iter = iter +1
		print 'invpower error = ',error
	x = x*conj(x[0])
	print 'Eval = ',e
	print 'Evect Real = ',norm(real(x))
	print 'Evect Imag = ',norm(imag(x))
	return x

def inversiter(eguess,wfguess,A,S):

    vk=np.matrix(wfguess)
    olu,opi=la.lu_factor(Ovr)
    
    for L in range(Klu):
        # get the LU decompostion of A-e*S
        lu,piv=la.lu_factor(A-eguess*S)
        
        for K in range(Kpower):
            vk=np.matrix(la.lu_solve((lu,piv),vk,overwrite_b=True))
            ek=(vk.T*A*vk).diagonal()/(vk.T*S*vk)
        
        if (eguess-ek[0,0])/eguess<1.e-9: return ek[0,0],wk[:,0]
        eguess=ek[0,0]
        vk=S*vk

def time_LU():
    """Print the times it takes to solve a system of equations using
    LU decomposition and (A^-1)B where A is 1000x1000 and B is 1000x500."""
    A = np.random.random((1000,1000))
    b = np.random.random((1000,500))
    start = time.time()
    L = la.lu_factor(A)
    a = time.time() - start
    start = time.time()
    A_inv = la.inv(A)
    a2 = time.time() - start
    start = time.time()
    la.lu_solve(L,b)
    a3 = time.time() - start
    start = time.time()
    np.dot(A_inv, b)
    a4 = time.time() - start

    
    time_lu_factor = a  # set this to the time it takes to perform la.lu_factor(A)
    time_inv = a2 # set this to the time it takes to take the inverse of A
    time_lu_solve = a3  # set this to the time it takes to perform la.lu_solve()
    time_inv_solve = a4  # set this to the time it take to perform (A^-1)B


    print "LU solve: " + str(time_lu_factor + time_lu_solve)
    print "Inv solve: " + str(time_inv + time_inv_solve)
    
    # What can you conclude about the more efficient way to solve linear systems?
    print "Better to use LU decomposition than inverse, cause it is NEVER a good idea to calculate an inerse"  # print your answer here.

 def updatedata(self, A):
     # Update b, c
     try:
         ALU = linalg.lu_factor(A)
         BC = linalg.lu_solve(ALU, c_[linalg.lu_solve(ALU, self.data.b + 1e-8*self.data.Brand[:,:1]), \
                              self.data.c + 1e-8*self.data.Crand[:,:1]], trans=1)
         C = linalg.lu_solve(ALU, BC[:,-1:])
         B = BC[:,:1]
     except:
         if self.C.verbosity >= 1:
             print 'Warning: Problem updating border vectors.  Using svd...'
         U, S, Vh = linalg.svd(A)
         B = U[:,-1:]
         C = num_transpose(Vh)[:,-1:]
 
     bmult = cmult = 1
     if matrixmultiply(transpose(self.data.b), B) < 0:
         bmult = -1
     if matrixmultiply(transpose(self.data.c), C) < 0:
         cmult = -1
     self.data.b = bmult*B*(linalg.norm(A,1)/linalg.norm(B))
     self.data.c = cmult*C*(linalg.norm(A,Inf)/linalg.norm(C))
     
     # Update
     if self.update:
         self.data.B[:,0] = self.data.b*(linalg.norm(A,1)/linalg.norm(self.data.b))
         self.data.C[:,0] = self.data.c*(linalg.norm(A,Inf)/linalg.norm(self.data.c))
         
         self.data.B[:,1] = self.data.w[:,2]*(linalg.norm(A,1)/linalg.norm(self.data.w,1))
         self.data.C[:,1] = self.data.v[:,2]*(linalg.norm(A,Inf)/linalg.norm(self.data.v,1))
         
         self.data.D[0,1] = self.data.g[0,1]
         self.data.D[1,0] = self.data.g[1,0]

    def linearize(self, inputs, outputs, partials):
        system_size = self.system_size
        self.lu = lu_factor(inputs['mtx'])

        partials['circulations', 'circulations'] = inputs['mtx'].flatten()
        partials['circulations', 'mtx'] = \
            np.outer(np.ones(system_size), outputs['circulations']).flatten()

def lnlike(h, X, y, covfunc):
    y = np.matrix(np.array(y).flatten()).T
    K = covfunc(X, X, h, wn = True)
    sign, logdetK = np.linalg.slogdet(K)
    alpha = np.mat(la.lu_solve(la.lu_factor(K),y))
    logL = -0.5*y.T * alpha - 0.5*logdetK - (y.size/2.)*np.log(2)
    return np.array(logL).flatten()

 def test_simple_known(self):
     # Ticket #1458
     for order in ['C', 'F']:
         A = np.array([[2, 1],[0, 1.]], order=order)
         LU, P = lu_factor(A)
         assert_array_almost_equal(LU, np.array([[2, 1], [0, 1]]))
         assert_array_equal(P, np.array([0, 1]))

def LU_det(A):
    (lu,piv) = la.lu_factor(A)
    
    # determine whether an even or odd number of row swaps
    s = (piv != np.arange(A.shape[0])).sum() % 2
    
    return ((-1)**s) * lu.diagonal().prod()

 def factored(self):
     "Caches the LU factorisation of the matrix"
     try:
         return self.lu_factored
     except AttributeError:
         self.lu_factored = la.lu_factor(self.val().simple_view())
         return self.lu_factored

def marglike(x, y, hyper, white_noise = False): # FIXME: build optional white noise into this kernel

    # Calculate covariance matrix
    K = matrifysquare(hyper, x, 0)
    K = 0.5* ( K + K.T) # Forces K to be perfectly symmetric

    # Calculate derivatives
    # dKdsigma = matrifysquare(hyper, x, 1) # Derivative w.r.t. log(sigma)
    # dKdlambda1 = matrifysquare(hyper, x, 2) # Derivative w.r.t. log(lambda1)
    # dKdh1 = matrifysquare(hyper, x, 3) # Derivative w.r.t. log(h1)

    sign, logdetK = np.linalg.slogdet(K)
    
    invKy = -0.5 * y.T *  np.mat(la.lu_solve(la.lu_factor(K),y)) \
        - 0.5 * logdetK - (y.size/2.) * np.log(2*np.pi)
    
    U = np.linalg.cholesky(K)

    n = len(x)
    L = - sum(np.log(np.diag(U))) -0.5 * y * invKy - n*0.5*np.log(2*np.pi)
    # dLdsigma = 0.5 * sum(np.diag(invKy*invKy.T*dKdsigma - (np.linalg.solve(K, dKdsigma)) ))
    # dLdlambda1 = 0.5 * sum(np.diag(invKy*invKy.T*dKdlambda1 - (np.linalg.solve(K, dKdlambda1)) ))
    # dKdh1 = 0.5 * sum(np.diag(invKy*invKy.T*dKdh1 - (np.linalg.solve(K, dKdh1)) ))

    return -L #, [-dKdsigma, -dKdlambda1, -dKdh1]

    def SolveNextTime(self):
        r"""
        Calculate the next time (factorization)
        and update the time stack grids
        """
        
        try:
            self.tstep += 1
        except :
            self.tstep = 0
            self.LinearSystem()
            # gets the m factor from the solved system
            self.mUtfactor = ln.lu_factor(self.mUt)

        self.Source(self.tstep)
        # in time
        # As t is in [0, 1, 2] (2nd order)
        # time t in this case is Utime[2]
        v = self.Independent()
        result = ln.lu_solve(self.mUtfactor, v)
        # reshape the vector to became a matrix again
        self.Ufuture = result
        # make the update in the time stack
        # before [t-2, t-1,  t]
        # after  [t-1, t,  t+1]
        # so t-2 receive t-1 and etc.
        self.Uprevious = self.Ucurrent 
        self.Ucurrent = self.Ufuture
        
        return self.Ufuture

def test_lapack_misaligned():
    M = np.eye(10,dtype=float)
    R = np.arange(100)
    R.shape = 10,10
    S = np.arange(20000,dtype=np.uint8)
    S = np.frombuffer(S.data, offset=4, count=100, dtype=np.float)
    S.shape = 10, 10
    b = np.ones(10)
    v = np.ones(3,dtype=float)
    LU, piv = lu_factor(S)
    for (func, args, kwargs) in [
            (eig,(S,),dict(overwrite_a=True)), # crash
            (eigvals,(S,),dict(overwrite_a=True)), # no crash
            (lu,(S,),dict(overwrite_a=True)), # no crash
            (lu_factor,(S,),dict(overwrite_a=True)), # no crash
            (lu_solve,((LU,piv),b),dict(overwrite_b=True)),
            (solve,(S,b),dict(overwrite_a=True,overwrite_b=True)),
            (svd,(M,),dict(overwrite_a=True)), # no crash
            (svd,(R,),dict(overwrite_a=True)), # no crash
            (svd,(S,),dict(overwrite_a=True)), # crash
            (svdvals,(S,),dict()), # no crash
            (svdvals,(S,),dict(overwrite_a=True)), #crash
            (cholesky,(M,),dict(overwrite_a=True)), # no crash
            (qr,(S,),dict(overwrite_a=True)), # crash
            (rq,(S,),dict(overwrite_a=True)), # crash
            (hessenberg,(S,),dict(overwrite_a=True)), # crash
            (schur,(S,),dict(overwrite_a=True)), # crash
            ]:
        yield check_lapack_misaligned, func, args, kwargs

    def SolveNextTime(self):
        r"""
        Calculate the next time (factorization)
        and update the time stack grids
        """

        try:
            self.tstep += 1
        except :
            self.tstep = 0
            self.LinearSystem()
            # gets the m factor from the solved system
            self.mUtfactor = ln.lu_factor(self.mUt)

        # As t is in [0, 1, 2] (2nd order)
        # time t in this case is Utime[2]
        # the independent term of the matrix, due the pressure field
        v = self.Independent()

        result = ln.lu_solve(self.mUtfactor, v)
        # reshape the vector to become a matrix again
        self.Ufuture = np.reshape(result, (self.Nz, self.Nx))

        # make the update in the time stack
        # before [t-2, t-1,  t]
        # after  [t-1, t,  t+1]
        # so t-2 receive t-1 and etc.
        # make the update in the time stack
        self.Uprevious[:][:] = self.Ucurrent[:][:]
        self.Ucurrent[:][:] = self.Ufuture[:][:]        
        
        return self.Ufuture