Have you tryed a simple newton algorithm (with the constraint added to the algo)? 


Let $(\alpha_{k})$ be defined as: $\alpha_k=1/k^2$  
  

> **Initialisation** : 
> 
> $x^0=[0,...,0]$
> 
> Compute $H=(A^* A)^{-1}$
> 
> **Loop for** $k$ in  $1:m$
> 
>   $x^k=x^{k-1}-\alpha_k  H A^*(Ax^{k-1}-b)$
> 
>   $x^{k}=g*x^{k}/\|x_2^k\|$
> 
> **end for loop**

Obviously, there are more adaptive way of choosing $\alpha_k$... but maybe you don't need such sofistication to solve a norm minimization problem. If $A^* A$ has very small eigen values you can use $H_k=(A^* A+\epsilon_k)^{-1}$ instead of $H$ ($\epsilon_k$ decreasing to zero)...


Note that this type of code is relatively general when you want to find the saddle point of a lagrangian and you know how to find the maxima with respect to Lagrange multipliers (in the dual space) (second step of the loop) but you need a gradient descent (or here Newton algo) for finding the minima in the principal space. 


Here is the corresponding R code: 


    A=t(array(1:1000,c(10,100)))
    m=100; b=1:10; g=3; l=5; p=10;
    alpha=1:m
    alpha=1/alpha^2
    x=array(0,c(m,p))
    H=t(A)%*%A
    svdH=svd(H)
    H=svdH$v%*%diag(1/svdH$d)%*%t(svdH$u) 
    for (k in 2:m)
    {
        x[k,]=x[k-1,]-alpha[k]*H%*%(t(A)%*%(A%*%x[k-1,]-b))
        x[k,]=g*x[k,]/sqrt(sum(x[k,(l+1):p]^2))
        print(sum((A%*%x[k,]-b)^2))
    }