Index: /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26685)
+++ /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26686)
@@ -138,5 +138,5 @@
 	return pe
 end #}}}
-function GetSolutionFromInputs(analysis::StressbalanceAnalysis,ug::Vector{Float64},element::Tria) #{{{
+function GetSolutionFromInputs(analysis::StressbalanceAnalysis,ug::IssmVector,element::Tria) #{{{
 
 	#Get dofs for this finite element
@@ -156,7 +156,7 @@
 
 		count += 1
-		ug[doflist[count]] = vx
+		ug.vector[doflist[count]] = vx
 		count += 1
-		ug[doflist[count]] = vy
+		ug.vector[doflist[count]] = vy
 	end
 
Index: /issm/trunk-jpl/src/jl/solve/modules.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/modules.jl	(revision 26685)
+++ /issm/trunk-jpl/src/jl/solve/modules.jl	(revision 26686)
@@ -125,5 +125,5 @@
 
 	#Initialize solution vector
-	ug = Vector{Float64}(undef,gsize)
+	ug = IssmVector(gsize)
 
 	#Go through elements and plug in solution
@@ -135,9 +135,9 @@
 
 end#}}}
-function InputUpdateFromSolutionx(analysis::Analysis,ug::Vector{Float64},femmodel::FemModel) #{{{
+function InputUpdateFromSolutionx(analysis::Analysis,ug::IssmVector,femmodel::FemModel) #{{{
 
 	#Go through elements and plug in solution
 	for i=1:length(femmodel.elements)
-		InputUpdateFromSolution(analysis,ug,femmodel.elements[i])
+		InputUpdateFromSolution(analysis,ug.vector,femmodel.elements[i])
 	end
 
@@ -145,5 +145,5 @@
 
 end#}}}
-function Reducevectorgtofx(ug::Vector{Float64},nodes::Vector{Node}) #{{{
+function Reducevectorgtofx(ug::IssmVector,nodes::Vector{Node}) #{{{
 
 	#Get size of output vector
@@ -151,13 +151,31 @@
 
 	#Initialize output vector
-	uf = Vector{Float64}(undef,fsize)
+	uf = IssmVector(fsize)
 
 	#Go through elements and plug in solution
 	for i=1:length(nodes)
-		VecReduce(nodes[i],ug,uf)
+		VecReduce(nodes[i],ug.vector,uf)
 	end
 
 	return ug
 
+end#}}}
+function Reduceloadx!(pf::IssmVector, Kfs::IssmMatrix, ys::IssmVector) #{{{
+
+	#Is there anything to do?
+	m, n = GetSize(Kfs)
+
+	if(m*n>0)
+
+		#Allocate Kfs*ys
+		Kfsy_s=IssmVector(m)
+
+		#Perform multiplication
+		MatMult!(Kfs,ys,Kfsy_s)
+
+		#Subtract Kfs*ys from pf
+		AXPY!(pf,-1.0,Kfsy_s)
+
+	end
 end#}}}
 function SystemMatricesx(femmodel::FemModel,analysis::Analysis)# {{{
Index: /issm/trunk-jpl/src/jl/solve/nodes.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/nodes.jl	(revision 26685)
+++ /issm/trunk-jpl/src/jl/solve/nodes.jl	(revision 26686)
@@ -136,9 +136,9 @@
 
 end# }}}
-function VecReduce(node::Node,ug::Vector{Float64},uf::Vector{Float64}) #{{{
+function VecReduce(node::Node,ug::Vector{Float64},uf::IssmVector) #{{{
 
 	for i=1:node.gsize
 		if node.fdoflist[i]!=-1
-			uf[node.fdoflist[i]] = ug[node.gdoflist[i]]
+			uf.vector[node.fdoflist[i]] = ug[node.gdoflist[i]]
 		end
 	end
Index: /issm/trunk-jpl/src/jl/solve/solutionsequences.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/solutionsequences.jl	(revision 26685)
+++ /issm/trunk-jpl/src/jl/solve/solutionsequences.jl	(revision 26686)
@@ -20,6 +20,12 @@
 		#Enforce constraints
 		ys = CreateNodalConstraintsx(femmodel.nodes)
+		Reduceloadx!(pf, Kfs, ys)
 
-		error("not implemented yet")
+		#Solve!
+		old_uf = uf
+		uf = Solverx(Kff, pf, old_uf)
+
+		print(uf)
+		error("compare with ISSM...")
 
 		#Increase count
Index: /issm/trunk-jpl/src/jl/solve/toolkits.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/toolkits.jl	(revision 26685)
+++ /issm/trunk-jpl/src/jl/solve/toolkits.jl	(revision 26686)
@@ -2,4 +2,6 @@
 #Toolkit #1: serial sparse arrays
 using SparseArrays
+
+#Matrix
 mutable struct IssmMatrix #{{{
 	matrix::SparseMatrixCSC{Float64,Int64}
@@ -20,5 +22,11 @@
 
 end#}}}
+function GetSize(matrix::IssmMatrix)#{{{
 
+	return size(matrix.matrix)
+
+end#}}}
+
+#Vector
 mutable struct IssmVector #{{{
 	vector::Vector{Float64}
@@ -26,4 +34,9 @@
 function IssmVector(M::Int64)#{{{
 	return IssmVector(zeros(M))
+end#}}}
+function GetSize(vector::IssmVector)#{{{
+
+	return length(vector.vector)
+
 end#}}}
 function AddValues!(vector::IssmVector,m::Int64,midx::Vector{Int64},values::Vector{Float64})#{{{
@@ -45,2 +58,24 @@
 
 end#}}}
+
+#Operations
+function MatMult!(A::IssmMatrix,x::IssmVector,y::IssmVector) #{{{
+
+	y.vector = A.matrix*x.vector
+
+end#}}}
+function AXPY!(y::IssmVector,alpha::Float64,x::IssmVector) #{{{
+
+	y.vector = alpha*x.vector + y.vector
+
+end#}}}
+function Solverx(A::IssmMatrix,b::IssmVector,xold::IssmVector) #{{{
+
+	x = IssmVector(GetSize(xold))
+	x.vector = (b.vector\A.matrix)' #not sure why we need to transpose... that's a mystery!
+
+	return x
+
+
+end#}}}
+