c
c
c***********************************************************************
c
      subroutine ORFAC1 ( kgp1, mmax, m, n, iprint, b, prjres, ierr )
c
c------------------------------------------------------------------------
c
c         stage 1:  householder reduction of
c
c                (    .    )      ( dr'. r3 )    kg
c                ( dr . r2 )  to  (----. -- ),
c                (    .    )      (  0 . r4 )  m-n-kg
c
c                  kg    1          kg   1
c
c      where dr = -d(q2)*y is the derivative of the modi-
c      fied residual produced by dpa, r2 is the trans-
c      formed residual from dpa, and dr' is in upper
c      triangular form (as in ref. (2), p. 18).
c      dr is stored in rows n+1 to m and columns n+2 to
c      n + kg + 1 of the matrix a (i.e.,columns 1 to kg
c      of the matrix b).  r2 is stored in column
c      n + kg + 2 of the matrix a (column kg + 1 of b).
c      for k = 1, 2, ..., kg, find reflection
c      i - u * u' / beta which zeroes b(i, k),
c      i = n+k+1, ..., m.
c
c  ......................................................
c
      real b(mmax,kgp1)
c
      kg=kgp1-1
      kg23=2*kg+3
      np1=n+1
      npkg=n+kg
c
      do 990 k=1,kg 
         alpha=0.0
         npk=n+k
         do 950 i=npk,m
950         alpha=alpha+b(i,k)**2
         alpha=sign(sqrt(alpha),b(npk,k))
         u1=b(npk,k)+alpha
         b(npk,k)=u1
         beta=alpha*u1
         if (alpha.ne.0.0) go to 960
c                          column was zero
         ierr=-8
c       ---------------------------------
         call VARERR (iprint,ierr,np1+k)
c       ---------------------------------
         go to 1040 
c                    apply reflections to remaining columns 
c                    of b and to residual vector. 
960      kp1=k+1
         do 980 j=kp1,kgp1
            acum=0.0
            do 970 i=npk,m
970            acum=acum+b(i,k)*b(i,j)
            acum=acum/beta
            do 980 i=npk,m
980            b(i,j)=b(i,j)-b(i,k)*acum
990      b(npk,k)=-alpha
c
      prjres=0.
      do 1000 j=np1,npkg
1000     prjres=prjres+b(j,kgp1)**2
      prjres=sqrt(prjres)
c
c        save upper triangular form and transformed
c        residual, for use in case a step is retracted.
c        also compute column lengths.
c
      if (ierr.eq.4) go to 1040
      do 1030 k=1,kg
         npk=n+k
         acum=0.
         do 1010 j=np1,npk
1010        acum=acum+b(j,k)**2
         do 1020 j=k,kgp1
            jsub=kgp1+j
            b(k,j)=b(npk,j)
1020        b(jsub,k)=b(npk,j)
         b(kg23,k)=sqrt(acum) 
1030  continue
c
1040  return
c
      end 
