solve.QP.compact <- function(Dmat, dvec, Amat, Aind, bvec, meq=0,
                   miter=1000, factorized=F){ 
  storage.mode(Dmat) <- "double"
  storage.mode(Amat) <- "double"
  storage.mode(Aind) <- "integer"
  storage.mode(dvec) <- "double"
  n     <- as.integer(nrow(Dmat))
  q     <- as.integer(ncol(Amat))
  anrow <- as.integer(nrow(Amat))
  if( missing(bvec) )
    bvec <- vector("double",q)
  else
    storage.mode(bvec) <- "double"
  iact  <- vector("integer",q)
  nact  <- 0
  r     <- min(n,q)
  sol   <- vector("double",n)
  crval <- 0
  work  <- vector("double",2*n+r*(r+5)/2+2*q+1)
  iter  <- vector("integer",2)

  res1 <- .Fortran("qpgen1", Dmat, dvec=dvec, n, n, sol=sol,
                   crval=as.double(crval),
                   Amat, Aind, bvec, anrow, q, as.integer(meq),
                   iact=iact, nact=as.integer(nact),
                   iter=as.integer(iter), 
                   as.integer(miter), work=work,
                   ierr=as.integer(factorized))

  if( res1$ierr == 1)
    stop("constraints are inconsistent, no solution!")
  else if( res1$ierr == 2)
    stop("maximal number of iteration exceeded.")
  else if( res1$ierr == 3)
    stop("matrix D in quadratic function is singular!")
    
  list(solution=res1$sol,
       value=res1$crval,
       unconstrained.solution=res1$dvec,
       iterations=res1$iter,
       iact=res1$iact[1:res1$nact])   
}

solve.QP <- function(Dmat, dvec, Amat, bvec, meq=0, 
                   miter=1000, factorized=F){ 
  storage.mode(Dmat) <- "double"
  storage.mode(Amat) <- "double"
  storage.mode(dvec) <- "double"
  n     <- as.integer(nrow(Dmat))
  q     <- as.integer(ncol(Amat))
  anrow <- as.integer(nrow(Amat))
  if( missing(bvec) )
    bvec <- vector("double",q)
  else
    storage.mode(bvec) <- "double"
  iact  <- vector("integer",q)
  nact  <- 0
  r     <- min(n,q)
  sol   <- vector("double",n)
  crval <- 0
  work  <- vector("double",2*n+r*(r+5)/2+2*q+1)
  iter  <- vector("integer",2)

  res1 <- .Fortran("qpgen2", Dmat, dvec=dvec, n, n, sol=sol,
                   crval=as.double(crval),
                   Amat, bvec, anrow, q, as.integer(meq),
                   iact=iact, nact=as.integer(nact),
                   iter=as.integer(iter), 
                   as.integer(miter), work=work,
                   ierr=as.integer(factorized))

  if( res1$ierr == 1)
    stop("constraints are inconsistent, no solution!")
  else if( res1$ierr == 2)
    stop("maximal number of iteration exceeded.")
  else if( res1$ierr == 3)
    stop("matrix D in quadratic function is singular!")

  list(solution=res1$sol,
       value=res1$crval,
       unconstrainted.solution=res1$dvec,
       iterations=res1$iter,
       iact=res1$iact[1:res1$nact])   
}



library.dynam("quadprog")
