$
$   rnormala.v2001   -   1/31/2001
$   rnormala.v2001   -   1/22/2001
$
$ file rnormala.v707 mag  last revised January 14, 2000
$ driver for rnormala subdmap, which processes SUPORT entries for the LAN 
$ method in a manner similar to that of the AHOU method.  Set up for SOLs 103,
$ 110-112
Compile phase1c nolist noref $
Alter 'type.*getlll' $ turn off L-set processing
getlll = false $ l-set not needed                               

Compile moders nolist noref $
Alter ' IF ( NORSET'(2,),''(2,) $ disable test for error 4407.  
$ Singular mrr allowed with alter.
If (false) then $ never enter this branch
Alter 'xread.*kxx' $ 
Type parm,,I,n,nr $ number of rigid body modes determined by alter
Type parm,,logical,n,rsetok $ returned false when number of rigid body
$ modes does not match size of r-set
Call rnormal phix,eigvmat,mxx,uset/pha1,phiac,phixx/s,nr/0.0/s,rsetok $ $
Equivx pha1/dra/always $ should be same as those computed by sekrrs
If (not rsetok) then $
  Message //'UFM.  number of rigid body modes is not equal to the'/
  ' number of DOFS listed on SUPORT entries'. $
  exit $
endif $
$ 'input vectors, modified outputs follow' $
Equivx phixx/phix/always $ replace original vectors with r-set refined vectors

$ file rna.v707
Compile rnormal nolist noref $
Subdmap rnormal phia,ruta,maa,uset/pha1,phiac,pha/nr/f2/rsetok $
$ Separate rigid body modes from phia into pha2, for rigid body modes, and 
$ phiac, for flexible modes. Uncouple the r-set modes of pha2 to make pha3, 
$ then reorthogonalize into pha1 such that pha1'*maa*pha1=i33.  Output the 
$ flexible modes phiac, then the appended modes pha=pha1|phiac.  
$ When there are more or less rb modes than the user prescribed, the input
$ modes are returned in pha, and pha1 and phiac represent the true rigid 
$ body modes and their complement, without r-set refinement.
$ All or any outputs may be purged.
Type parm,,I,n,norset,nr1,nc,c126,nr $ nr is number of r.b. modes, returned as 
$ -1 when there appear to be none. output
Type parm,,rs,n,fi,f2,rut2 $ f2 is a root larger than the expected size of a  $ zero root.  If input as 0.0 it is reset to larger number.  Optional input.
$ In practice, set to f2 on eigrl, if defined.
Type parm,,I,n,ii $ loop counter
Type parm,,rs,y,rnratio=1.e6 $ ratio between f2 and largest freq. accepted
$ as almost-zero root.
Type parm,,logical,n,rsetok $ returned true when nr=norset.  Output
Paraml uset//'uset'//////'r'/s,n,norset $
Paraml phia//'trailer'/1/s,n,nr1 $ number of all modes computed
If (f2=0.0) then $
  Norm ruta/scra///s,n,rut2 $ largest root
  If (rut2<0.1) rut2=0.1 $ guard against lots of rb modes
Else $
  Rut2 = (pi(2)*f2)**2
Endif $
Rut2 = rut2/rnratio $
Do while (ii<nr1) $
  Ii = ii + 1 $
  Paraml ruta//'dmi'/ii/ii/s,n,ruti
  If (ruti>rut2) then $
    Ii = ii-1 $
    Jump loopend $
  Endif $
Enddo $
Label loopend $
Nr = ii $ $ reset
Fi = sqrt(ruti)/pi(2) $
Message //'uim.  first flexible mode frequency is '/fi/'hz.' $
Message //'uim.  number of rigid body modes is'/nr $   
IF (NR><NORSET) JUMP BADRSET $         
Nc = nr1-nr $
Matgen ,/v1r/6/nr1/nr/nc $ 
if (nc>0) then $
  Partn phia,v1r,/pha2,,phiacx,/0 $ arbitrary rigid body modes
  Equivx phiacx/phiac/always $ allow purged output
Else $
  Equivx phia/pha2/always $
Endif $
if (norset>-1 and nr=norset) then $
  upartn uset,pha2/phrr,,,/'a'/'r'/'l'/1 $ 
  CEAD  phrr,,,,,,/vt,CLAMA,OCEIGS,u,S/
    S,N,NFOUND//-1/'SVD'//-1000 $ FIND up to 1000 SINGULAR VALUES
  $ phrr = u*s*vt', phrr-1 = vt*sinv*u'
  Paraml S//'dmi'/1/1/s,n,maxrut $ roots are gen. mass
  Paraml S//'dmi'/nr/nr/s,n,minrut $ ordered largest first
  If (minrut*rnratio < maxrut) jump badrset $
  solve s,,,,/sinv/3 $ 
  trnsp u/ut
  smpyad vt,sinv,ut,,,/phinvrr/3 $
  $ofp clama// $
  mpyad pha2,phinvrr,/pha3 $
  Mpyad maa,pha3,/mpha3 $ 
  Mpyad pha3,mpha3,/x22/1////6 $ gen. Mass in arbitrary basis
  Trnsp pha3/pha3t $
  C126=GETSYS(C126,126) $ CAPTURE SYS(126)  method to turn off sparse methods                                
  PUTSYS(0,126) $                                                           
  Decomp x22/l22,,//1 $ 1= CHOLESKY, X=L'*L
  Fbs l22,,pha3t/pha1t///1 $ 1=fwd only
  PUTSYS(C126,126) $  reset to orignal value
  Trnsp pha1t/pha1 $
  $Smpyad pha1,maa,pha1,,,/genmchk/3////1 $
  $Matprn genmchk// $
  Message //'uim.  Modes have been refined by r-set operations' $
Else $
  Equivx pha2/pha1/always $ put number of rbmodes in first output,
  $ based on size of # zero roots, not size of r-set
  Jump badrset $
Endif $
jump lok $
Label badrset $ note that there are two ways to get here
Equivx phia/pha/always $ equiv input to output
If (norset><nr) message //'uwm.  Size of r-set is'/norset/
' Number of rigid body modes is'/nr $
message //'Mode refinement by r-set branched over.'/
 'adjust r-set (suport entries) to match zero. Freq. modes' $ $
rsetok = false $
return $
label lok $
if (nc>0) then $
  merge  pha1,,phiacx,,v1r,/phax/0 $ refined r.b. modes, flex modes
  equivx phax/pha/always $ allow purged output
Else $
  Equivx pha1/pha/always $
Endif $
Rsetok = true $
Return $ 
End $ rnormal