ginv.inc

心理测量统计分析中

(*-----------------------------------------------------

            Generalized Inverse

    When
         RAC = I-              -I = J
               I  1 0 ...     0 I
               I  0 1 0 ...   0 I
               I      .         I
               I      .         I
               I      .         I
               I      0 1 0 ... I
               I      0 0 0 ... I
               I      .         I
               I      .         I
               I      .         I
               I-              -I
    then
         g-Inverse of A = CJ'R

-------------------------------------------------------*)

//type  TMatGInv = array[1..20,1..20] of Extended;

procedure GInv( a : TMatGInv;       //   m * n matrix
                var b : TMatGInv;   //   n * m matrix
                m, n  : Longint;
                var rank : Longint );
  Label QP;
  const acc = 1.0e-15;
  var r, c : TMatGInv;
      i, j, i1, j1, ki, kj : Longint;
      CkTranspose : (Yes, No);
      zero, v : Extended;

  procedure Swap( var a, b : Extended );
    var t : Extended;
    begin
          t:=a;  a:=b;  b:=t;
    end;

  begin   {   GInv   }
      if m <= n
        then CkTranspose:=No
        else
          begin
            CkTranspose:=Yes;
            r:=a;
            for i:=1 to n do
              for j:=1 to m do
                a[i,j]:=r[j,i];
            i:=m; m:=n; n:=i;
          end;

      v:=0.0;
      for i:=1 to m do
        for j:=1 to n do
          if v < abs(a[i,j]) then v:=abs(a[i,j]);
      if v > acc
        then zero:=acc*v
        else
          begin
              rank:=0;
              goto QP;
          end;

      for i:=1 to m do
        for j:=1 to m do
          if i = j then r[i,j]:=1.0
                   else r[i,j]:=0.0;
      for i:=1 to n do
        for j:=1 to n do
          if i = j then c[i,j]:=1.0
                   else c[i,j]:=0.0;

      rank:=m;
      for i:=1 to m do
        begin
          ki:=i; kj:=i;
          v:=abs(a[ki,kj]);
          for i1:=i to m do
            for j1:=i to n do
              if v < abs(a[i1,j1]) then
                begin
                  ki:=i1; kj:=j1;
                  v:=abs(a[ki,kj]);
                end;
          if v < zero
            then
              begin
                rank:=i-1;
                goto QP;
              end
            else
              begin
                if ki > i then
                  begin
                    for j1:=i to n do
                      Swap( a[i,j1], a[ki,j1] );
                    for j1:=1 to m do
                      Swap( r[i,j1], r[ki,j1] );
                  end;

                if kj > i then
                  begin
                    for j1:=i to m do
                      Swap( a[j1,i], a[j1,kj] );
                    for j1:=1 to n do
                      Swap( c[j1,i], c[j1,kj] );
                  end;

                v:=1/a[i,i];
                for j1:=i to n do
                  a[i,j1]:=v*a[i,j1];
                for j1:=1 to m do
                  r[i,j1]:=v*r[i,j1];
                if i < m then
                  begin
                    for i1:=i+1 to m do
                      begin
                        v:=a[i1,i];
                        for j1:=i to n do
                          a[i1,j1]:=a[i1,j1]-v*a[i,j1];
                        for j1:=1 to m do
                          r[i1,j1]:=r[i1,j1]-v*r[i,j1];
                      end;
                  end;
                if i < n then
                  begin
                    for j1:=i+1 to n do
                      begin
                        v:=a[i,j1];
                        for i1:=i to m do
                          a[i1,j1]:=a[i1,j1]-v*a[i1,i];
                        for i1:=1 to n do
                          c[i1,j1]:=c[i1,j1]-v*c[i1,i];
                      end;
                  end;
              end;
        end;

    QP : ;
      for i:=1 to n do
        for j:=1 to m do
          begin
            v:=0.0;
            if rank > 0 then
              for i1:=1 to rank do
                v:=v+c[i,i1]*r[i1,j];    
            b[i,j]:=v;
          end;
      if CkTranspose = Yes then
        begin
          r:=b;
          i:=m; m:=n; n:=i;
          for i:=1 to n do
            for j:=1 to m do
              b[i,j]:=r[j,i];
        end;
  end;      {   GInv   }