$ nlirma.v2001  R. L. Harder, modified by mag 
$ mag changes in lower case.  Last revised June 13, 2003
$ options added for: 
$ Getting rigid body modes in eigensolution.
$ Moving the reference coordinates
$ Output of constraint forces as mpcforces
$ ALTER FOR INERTIA RELIEF IN SOL 106   BY BOB HARDER   05/11/03
$   NEEDS GRID AND USER SET U2 FOR LAGRANGE DOFS
$   DO NOT SPC RIGID BODY MODES (LAGRANGE MULTIPLIER CONSTRAINTS)
$   PARAM,TESTNEG IS SET TO 1 (CONTINUE IF NOT POSITIVE DEFINITE)
$
$ revision 1.  rgnl error fixed for mce1 input, failed experiment removed from 
$ comments.
$ECHOOFF
COMPILE NLSTATIC list $
Alter 'file.*uset1' $
FILE     rgnl=SAVE,OVRWRT $ now used in loop, save for next loop.
ALTER 'IF.*NEWK.*MSCHG','' $ ALWAYS: STIFFNESS GENERATION AND REDUCTION
IF ( MSCHG > -100 ) THEN $
ALTER 'IF.*NEWP.*MSCHG','' $ NEVER: SKIP LINEAR IF STILL SAME LOAD
IF ( MSCHG > -100 ) THEN $
$=====================================================================
ALTER 'DBVIEW *KGG0','' $
Paraml uset1//'uset'//////'u2'/s,n,nou2set/'u3'/s,n,nou3set $
VEC	USET1/VG12/'G'/'COMP'/'U2' $ LAGRANGE MULTIPLIERS
PARAML	VG12//'TRAILER'/5/S,N,L2SET $
PARAM	//'PREC'/S,N,PREC $
PARAM	//'DIV'/S,N,L2SET/V,N,L2SET/V,N,PREC $
IF ( L2SET <> 6 ) THEN
  MESSAGE //'ufm. SIX LAGRANGE DEGREES OF FREEDOM NEEDED IN USET U2' $
  EXIT $
ENDIF
$MESSAGE //' INERTIA RELIEF, SIX LAGRANGE MULTIPLIERS' $ prints on every loop
VECPLOT MUGNI,BGPDTS,,CSTMS,CASES,,,,/UGVBASU,,,,/0/0/3 $
MATMOD  UGVBASU,BGPDTS,,,,/BGPDTU,/11 $
VECPLOT,,BGPDTU,,CSTMS,,,,,/RU2G,,,,/grdpnt/0/4 $ UPDATED RIGID VECTORS
DBVIEW  KGGV=KGG (WHERE NLOOP=0) $
PARTN   KGGV,VG12,/K11,,, $
PARTN   MJJ ,VG12,/M11,,, $
PARTN   RU2G,VG12,/RU21,,,/1 $
MPYAD   RU21,M11,/L21 $
TRNSP	L21/L12 $
Trnsp   l21/l21t $
Mpyad   ru21,l21t,/m22/////6 $  
Call dbstore l21,m22,,,//0/0/'     '/0 $ save for before statics super3
MERGE   K11,L21,L12,,VG12,/KGG0 $
EQUIVX	KGG0/KGG0X/AlWAYS $ FOR V2004
$=====================================================================
alter 'dcmp.*kllrh'(,-1) $ get rid of terms on lmt diags
SETVAL  //S,N,TESTNEGX/1 $
VEC	USETNL/VL12/'L'/'COMP'/'U2' $ LAGRANGE MULTIPLIERS
PARTN   KLLRH,VL12,/K11RH,K21RH,K12RH, $
DELETE	/KLLRH,,,, $
MERGE	K11RH,K21RH,K12RH,,VL12,/KLLRH $
$=====================================================================
ALTER 'NLITER.*KAAL'(,-1) $
PARTN   KAAL,VL12,/K11L,K21L,K12L, $
DELETE	/KAAL,,,, $
MERGE	K11L,K21L,K12L,,VL12,/KAAL

Alter 'EQUIVX   KLLRH/KAA/-1','' $
Type parm,,char8,y,rbmodes='no' $
If (rbmodes='yes') then $
  matgen ,/i22/1/nou2set $
  VEC	USETNL/VL12x/'L'/'COMP'/'U2' $ LAGRANGE MULTIPLIERS
  PARTN   KLLRH,VL12x,/K11RH1,,, $
  MERGE	K11RH1,,,i22,VL12x,/kaa $
  $ discard the lmt coefficients
  Message //'uim.  Lmt coefficients removed before eigensolution' $
Else $
  EQUIVX   KLLRH/KAA/-1 $
Endif $

Alter 'super3.*uln'(,-1) $ just before data recovery.  Ugn is an input/output
$ remove the lmt variables from the output vector ugn
upartn usetnl,ugn/mq2,u1nx,,/'g'/'u2'//1 $ get the lmt solution for force
umerge1 usetnl,,u1nx,,/ugntemp/'g'/'u2'//1 $ merge in zero for lmt values
equivx ugntemp/ugn/always $ contains only displacements now

$ shift reference axes, if requested
Type parm,,I,y,scdir $ direction of scaling variable
Type parm,,I,n,nc,scdirm1 $ local variables
If (scdir<0 or scdir>3) then $
  Message //'ufm.  Param, scdir has a value of'/scdir $
  Message //'it must be in the range 0 to 3' $
  Exit $
Endif $
If (scdir>0) then $
  Paraml cstms//'presence'////s,n,nocstms $
  If (nocstms>-1) then 
    Message //'ufm.  the scdir option can not be used when cordi '/
      'entries are present.'
    Exit $
  Endif $
  Message //'uim. Selected reference direction is '/scdir $ 
  Matmod eqexins,usetnl,sils,casecc,,/vg1c,/17//1 $ has 1.0 for 6 dof partn
  $ of the GID listed for the PARTN case control command
  Nc = 6-scdir $
  Scdirm1 = scdir-1 $
  Matgen ,/u16/6/6/scdirm1/1/nc $ has 1.0 where selected
  Merge, ,u16,,,,vg1c/vg1d/0 $ has 1.0 at selected g-set location
  Message //'GID of selected reference node'
  matgpr gpls,usetnl,sils,vg1d//'h'/'g' $
  Partn ugn,,vg1d/,u1n,,/0 $ value at the selected DOF
  Paraml u1n//'trailer'/1/s,n,ncols $
  If (ncols>1) then $
    Trnsp u1n/u1nt $
    Matmod u1nt,,,,,/u11n,/28 $ diagonalize
    Matgen, /unitc/6/ncols/0/ncols $
    Trnsp unitc/unitct $
    Mpyad u16,unitct,/uc6 $ put same thing in each column
  Else $
    Equivx u1n/u11n/always $
    Equivx u16/uc6/always $
  Endif $
  Mpyad uc6,u11n,/uc6n $ all six values the same
  Matmod uc6n,sils,,,,/u1gn,/4 $ has the same value at all gids
  Upartn usetnl,u1gn/scrn,u12n,,/'g'/'u2'//1 $
  Umerge1 usetnl,,u12n,,/u2gn/'g'/'u2'//1 $ discard correction on lmt variable 
  add ugn,u2gn/ugntemp1//-1.0 $ subtract off the selected value everywhere
  equivx ugntemp1/ugn/always $
  matprn uc6n// $ label matrices after print to avoid extra page eject
  message //'uim.  The deflections of uc6n above are subtracted '/
    'from all dofs, at user request with param,scdir1' $
endif $

$ output the forces on the masses due to the constraints as mpc forces
call dbfetch /l21,m22, ,,/0/0/0/0/0 $ phit*m
mpyad l21,mq2,/qm1n/1/-1 $ m22*phi*(-q1) = mpcforces
umerge1 usetnl,,qm1n,,/qmgn1/'g'/'u2'//1 $
mpyad m22,mq2,/f2rb $
matprt m22// $
message //'rigid body mass matrix, above'
matprt mq2// $ because the constraints are mass weighted, mq2 is actually
$ an acceleration matrix, dimensionally.  It equals the resultants of the 
$ mpc forces below.
message //'rigid body acceleration vector, above' $
matprt f2rb// $ 
message //'generalized force on LMT variables, above' $
if (nomset=-1) then $ use mpcforce output slot
  message //'inertial forces are output as MPCFORCES for'/
  ' DOFs selected by mpcforce =' $
  call xadd5 qmgn1,qmgn,,,/qmgn2/cu/cu/cz/cz/cz $
  equivx qmgn2/qmgn/always $
else $ print with ofp-like vector printer
  matgpr bgpdts,usetnl,,qmgn1//'ofp'/'g'//////
    'load num'/'BER'///////                                                        
    'I n e r '/'t i a l '/' f o r c'/'e       ' $ 
endif $ end of nlirma alter package