Attribute VB_Name = "Harris_Kaiser"
Public nv, m_fact
Public Sub 回転Harris_Kaiser()
' power 0から1より小の値. 0の時独立解（相関が高すぎる） 0.5の時proportional solution.
   pw = 0  '1 / 3 '0.5 ' 1 / 3
Dim st$, j
For lno = 1 To 2
   If lno = 2 Then pw = 0.5
   j = 1
   st$ = " Factor pattern matrix, Communality and Unique variance"
'st$ = " Orthomax rotated factor pattern matrix, Communality and Unique variance"
   Do
     j = j + 1
   Loop Until (Cells(j, 1) = st$) Or (j > 10000)
   If j > 10000 Then
     MsgBox (st$ + "の解を見つけることができませんでした。終了します。")
     End
   End If
   j = j + 1 '+ 1

'因子数を調べる
   k = 2
   Do
     k = k + 1
   Loop Until (Cells(j, k) = "" Or k > 25)
   If k > 25 Then
     MsgBox ("因子数がわかりませんでした。終了します。")
     End
   End If
   m_fact = Cells(j, k - 1)
'   MsgBox m_Fact

'変数の数を数える
   k = j
   Do
     k = k + 1
   Loop Until (Cells(k, 1) = "Variance" Or Cells(k, 1) = "" Or k - j > 200)
   If (k - j) > 200 Then
     MsgBox ("因子数がわかりませんでした。終了します。")
     End
   End If
  nv = k - j - 1

  ReDim a_fact(nv, m_fact), varnum(nv), varnames(nv) As String
   For jj = 1 To nv
     For ii = 1 To m_fact
       a_fact(jj, ii) = Cells(j + jj, ii + 2)
     Next ii
     varnum(jj) = Cells(jj + j, 1)
     varnames(jj) = Cells(jj + j, 2)
   Next jj
' 書き込み開始セルを見つける。、

   istart = 10000
   Do
     istart = istart - 1
   Loop Until (Cells(istart, 1) <> "" Or istart <= 1)
   If istart <= 1 Then
     MsgBox "書き込み開始位置を決めることができませんでした。"
     End
   End If
   
   istart = istart + 1
   For jj = 1 To nv
     Cells(jj + istart + 2, 1) = varnum(jj)
     Cells(jj + istart + 2, 2) = varnames(jj)
   Next jj
   Cells(istart, 1).Activate
   Call HKOB(nv, m_fact, a_fact, 1, pw)
'   Call HKOB
   Call sorting(istart + 3, 3, m_fact + 3)
   Next lno
   End

End Sub

'  M=変数数, MF=因子数, D=負荷量、IR=設定している変数の数,
'  PW=Harris-Kaiser の係数　、IC=因子数（設定数）
Private Sub HKOB(m, mf, D(), INRM, pw)
'Sub HKOB()

'C
'C     HARRIS-KAISER CASE II ORTHOBLIQUE ROTATION WITH ASSOCIATED RESULTS
'C
'C     MODEL: COVARIANCE STRUCTURE IN COMMON FACTOR SPACE
'C                 =Q*EM*Q'
'C                 =Q*EM**((2-PW)/2)*T3*D'
'C
'C     Q (EIGEN VECTORS)         --- CONSTRAINED　　 固有ベクトル
'C     EM (EIGEN VALUES)         --- FREE PARAMETERS 固有値
'C     T3 (ROTATION MATRIX)       --- CONSTRAINED　　変換行列
'C     D (=Q*EM**(PW/2)*T3 : RESCALED ROTATED-LOADINGS) --- CONSTRAINED
'C
'C     D4 (ROTATED LOADINGS　回転後負荷量), PH (FACTOR CORRELATIONS　因子間相関),
'C     DI: SQRT OF DIAG( T3*EM**((2-PW)/2)*T3 )
'C     S1: CONTRIBUTIONS OF ORTHOMAX-ROTATED OBLIQUE-FACTORS BEFORE SCALING
'C     S1A: CONTRIBUTIONS OF ORTHOBLIQUE-ROTATED FACTORS　因子の寄与
'C     HS1: COMMUNALITIES 共通性
'C
'C     W:ORTHOMAX WEIGHT
'C     INRM---0:RAW ORTHOMAX, 1:NORMALIZED ORTHOMAX
'C        <HKOB>                         '99.10.8
'C
'      IMPLICIT REAL*8(A-H,O-Z)
IR = m
IC = mf
      ReDim HS1(IR), PC(IR, IR), q(IR, IR), EM(IR), U(IR, IC)
      ReDim ND(IR, IC), T31(IC, IC), T32(IC, IC), T3(IC, IC), PH(IC, IC), DI(IC)
      ReDim d3(IR, IC), d4(IR, IC), S1(IC), S1A(IC)
'C
'C     FACTOR ROTATION
      For i = 1 To m
        HS1(i) = 0
        For j = 1 To mf
          d3(i, j) = D(i, j)
          HS1(i) = HS1(i) + D(i, j) ^ 2
        Next j
      Next i
'      write(6,345)(hs1(i),i=1 to m)  共通性を書き出す
'  345 format(/,'COM=',12f7.3)
'C
      For i = 1 To m
        For j = 1 To i
          PC(i, j) = 0
          For k = 1 To mf
            PC(i, j) = PC(i, j) + D(i, k) * D(j, k)
            PC(j, i) = PC(i, j)
          Next k
        Next j
      Next i
      eps = 0.000001
      MAXI = 500
'C
'      Call EIGENx(PC, q, m, eps, MAXI, IR)
      Call JACOBI(PC, m, q, eps)
'C
      For i = 1 To mf
        EM(i) = PC(i, i)
      Next i
      For i = 1 To m
        For j = 1 To mf
          U(i, j) = q(i, j) * EM(j) ^ (pw / 2#)
          D(i, j) = U(i, j)
        Next j
      Next i
'C
      Call OTMAX(m, mf, D, S1, w, HS1, INRM, IR, IC)
'C
'C     ROTATION MATRIX
      For i = 1 To mf
        For j = 1 To mf
          T31(i, j) = 0
          T32(i, j) = 0
           For k = 1 To m
             T32(i, j) = T32(i, j) + U(k, i) * D(k, j)
             T31(i, j) = T31(i, j) + U(k, i) * U(k, j)
           Next k
        Next j
      Next i
'C
      Call INV(T31, mf, IC, DET)
'C
      For i = 1 To mf
        For j = 1 To mf
          T3(i, j) = 0
          For k = 1 To mf
            T3(i, j) = T3(i, j) + T31(i, k) * T32(k, j)
          Next k
        Next j
      Next i
'C
'C     ROTATED LOADINGS AND FACTOR CORRELATIONS
      For i = 1 To mf
        For j = 1 To mf
          PH(i, j) = 0
          For k = 1 To mf
            PH(i, j) = PH(i, j) + T3(k, i) * EM(k) ^ (1# - pw) * T3(k, j)
          Next k
        Next j
        DI(i) = Sqr(PH(i, i))
      Next i
      For i = 1 To mf
        For j = 1 To mf
          PH(i, j) = PH(i, j) / (DI(i) * DI(j))
        Next j
      Next i
      For i = 1 To m
        For j = 1 To mf
          d4(i, j) = D(i, j) * DI(j)
        Next j
      Next i
      For j = 1 To mf
        S1(j) = 0
        S1A(j) = 0
        For i = 1 To m
          S1(j) = S1(j) + D(i, j) ^ 2
          S1A(j) = S1A(j) + d4(i, j) ^ 2
        Next i
      Next j
'C
      ii = ActiveCell.Row
      ii = ii + 1
      Cells(ii, 1) = "ROTATED LOADINGS (HARRIS-KAISER ORTHOBLIQUE METHOD) Power= " + CStr(pw)
      Cells(ii, 1).Activate
      Call cellname(CInt(ii), 1, "_Harris_Kaiser因子パタン")
      ii = ii + 1
      For j = 1 To mf
         Cells(ii, j + 2) = j
      Next j
      For i = 1 To m
         For j = 1 To mf
           Cells(ii + i, j + 2) = d4(i, j)
            Cells(ii + i, j + 2).NumberFormatLocal = "###0.000_ "
            If Abs(Cells(ii + i, j + 2)) >= 0.4 Then
              Call coloring(ii + i, j + 2, 34)
             End If

'        WRITE(6,136)I,(D4(I,J),J=1 to MF)
         Next j
'  136 FORMAT(I6,2X,7F8.4)
      Next i
    Range(Cells(ii + m, 1), Cells(ii + m, mf + 1 + 1)).Select
    With Selection.Borders(xlEdgeBottom)
        .LineStyle = xlDouble
        .Weight = xlThick
        .ColorIndex = xlAutomatic
    End With
      ii = ii + i
       Cells(ii, 1).Activate
      Cells(ii, 2) = "CONTR."
      For i = 1 To mf
         Cells(ii, i + 2) = S1A(i)
         Cells(ii, i + 2).NumberFormatLocal = "##0.000_ "
      
      Next i
'      WRITE(6,137)(S1A(I),I=1 to MF)
'  137 FORMAT(/'CONTR.=',1X,7F8.4)
      ii = ii + 2
      Cells(ii, 1) = "因子間相関(Harris-Kaiser)"
      Call cellname(CInt(ii), 1, "_Harris_Kaiser因子間相関")

      ii = ii + 1
      For i = 1 To mf
         Cells(ii, i + 2) = i
      Next i
      For i = 1 To mf
        Cells(ii + i, 2) = i
        Cells(ii + i, 1) = " "
        For j = 1 To mf
         Cells(ii + i, j + 2) = PH(i, j)
         'Cells(ii + i, j + 2).NumberFormatLocal = "###0.000_ "
         Call corrcolor(ii + i, j + 2)
        Next j
        ' WRITE(6,139)(PH(I,J),J=1 to I)
'         139 FORMAT(/'FAC.COR.=',7F8.4)
      Next i
 '     Return
      End Sub

'C
Private Sub OTMAX(m, mf, a, CT, w, CM, INRM, IR, IC)

'C
'C     ORTHOMAX ROTATION
'C     M=# OF VARIABLES, MF=# OF FACTORS,
'C     A:LOADING MATRIX, CT:CONTRIBUTIONS, CM:COMMUNALITIES
'C     W:ORTHOMAX WEIGHT
'C     INRM---0:RAW ORTHOMAX, 1:NORMALIZED ORTHOMAX
'C
'C     <OTMAX>                                            '99.9.3
'C
'      IMPLICIT REAL*8(A-H,O-Z)
   '   Dim A(IR, IC), CM(IR), CT(IC)
      PAI = 3.14159265358979
      eps = 0.000001
      imax = 500
      For i = 1 To m
        CM(i) = 0
        For j = 1 To mf
          CM(i) = CM(i) + a(i, j) ^ 2
        Next
      Next
'C
      If (INRM = 1) Then
        For i = 1 To m
          For j = 1 To mf
            a(i, j) = a(i, j) / Sqr(CM(i))
          Next j
        Next i
      End If
'C
      am = m
      it = 0
1     it = it + 1
      If (it > imax) Then
        MsgBox "otmax: NO CONVERGENCE : # OF ITERATIONS=" + Str(it - 1)
        End
      End If
'C
      IND = 0
      For i = 1 To mf - 1
        For j = i + 1 To mf
          aa = 0
          bb = 0
          cc = 0
          dd = 0
          For k = 1 To m
            E = a(k, i) ^ 2 - a(k, j) ^ 2
            F = 2# * a(k, i) * a(k, j)
            aa = aa + E
            bb = bb + F
            cc = cc + E ^ 2 - F ^ 2
            dd = dd + 2# * E * F
          Next k
'C
          t1 = dd - 2# * aa * bb / am * w
          t2 = cc - (aa ^ 2 - bb ^ 2) / am * w
          If (t2 = 0) Then
             MsgBox "TANGENT IS INFINITE !!!!"
             End
           End If
          T = Atn(t1 / t2)
          If (t2 < 0 And t1 >= 0) Then T = T + PAI
          If (t2 < 0 And t1 < 0) Then T = T - PAI
          T = T / 4
'C
          If (Abs(T) >= eps) Then IND = 1
          If (Abs(T) >= eps) Then
             x = Cos(T)
             y = Sin(T)
             For k = 1 To m
               z = a(k, i) * x + a(k, j) * y
               a(k, j) = -a(k, i) * y + a(k, j) * x
               a(k, i) = z
             Next k
          End If
        Next j
      Next i
      If (IND = 1) Then GoTo 1
'C
      For j = 1 To mf
        CT(j) = 0
        For i = 1 To m
          If (INRM = 1) Then a(i, j) = a(i, j) * Sqr(CM(i))
          CT(j) = CT(j) + a(i, j) ^ 2
        Next i
      Next j
'      WRITE(6,101)INRM,W,IT,EPS
'  101 FORMAT(//'ORTHOMAX SOLUTION: KAISER''S NORMALIZATION=',I1,
'     & ' (0:NO, 1:YES)',//,'ORHTOMAX WEIGHT=',F5.2//'# OF ITERATIONS=',
'     & I3,' EPS=',F11.7/)
'C
'      for I=1 to M
'        WRITE(6,103)I,(A(I,J),J=1 to MF)
'  103 FORMAT(I5,2X,(8F9.3))
'      next
'      WRITE(6,104)(CT(I),I=1 to MF)
'  104 FORMAT(/'CTRBTN=',(8F9.3))
 '     Return
 End Sub

      Sub EIGENx(b, w, m, E, MAXI, L)
'C     B:固有値（対角）, W:固有ベクトル, M:変数の数(ここで実際に処理する）
'           , E:精度 , MAXI:最大反復数 , L変数の数 （呼び出しもとの大きさ）
'C     EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC MATRIX
'C     B -> W B W', W'W=I, B: DIAGONAL, B(I,I)>=B(I+1,I+1), (I=1,M-1)
'C      <EIGEN>                                '98.4 '99.10.5
'C
'      IMPLICIT REAL*8(A-H,O-Z)
      ReDim a(L, L)    ',B(L,L),W(L,L)
      For i = 1 To m
        For j = 1 To m
          w(i, j) = 0
          If (i = j) Then w(i, j) = 1
        Next j
      Next i
'C
      ITER = 0
    Do
      c = Abs(b(2, 1))
      i1 = 2
      j1 = 1
      For i = 2 To m
        For j = 1 To i - 1
         If (Abs(b(i, j)) > c) Then
           i1 = i
           j1 = j
           c = Abs(b(i, j))
         End If
        Next j
      Next i
      If (c < E) Then Exit Do 'go to 5   '''''
      T = 2# * b(i1, j1) / (b(i1, i1) - b(j1, j1) + Sqr((b(i1, i1) - b(j1, j1)) ^ 2 + 4# * b(i1, j1) ^ 2) _
         * Sgn(b(i1, i1) - b(j1, j1)))    '---------------------------------
      p = 1# / Sqr(1# + T ^ 2)
      q = T * p
'C
      For i = 1 To m
        If (i <> i1 And i <> j1) Then
          x = b(i, i1) * p + b(i, j1) * q
          y = -b(i, i1) * q + b(i, j1) * p
          b(i, i1) = x
          b(i, j1) = y
          b(i1, i) = x
          b(j1, i) = y
        End If
        x = w(i, i1) * p + w(i, j1) * q
        y = -w(i, i1) * q + w(i, j1) * p
        w(i, i1) = x
        w(i, j1) = y
      Next i
'C
      BII = b(i1, i1)
      BJJ = b(j1, j1)
      BIJ = b(i1, j1)
      b(i1, i1) = BII * p ^ 2 + BJJ * q ^ 2 + 2# * BIJ * p * q
      b(j1, j1) = BII * q ^ 2 + BJJ * p ^ 2 - 2# * BIJ * p * q
      b(i1, j1) = (BJJ - BII) * p * q + BIJ * (p ^ 2 - q ^ 2)
      b(j1, i1) = b(i1, j1)
'C
      ITER = ITER + 1
      If (ITER > MAXI) Then
         MsgBox "固有値：NO CONVERGENCE !!! : NO. OF ITRNS= " + Str(ITER)
         End
      End If
'C
    Loop
      For i = 1 To m - 1
        i1 = i
        c = b(i, i)
        For j = i + 1 To m
          If (b(j, j) > c) Then
            i1 = j
            c = b(j, j)
          End If
        Next j
        If (i1 <> i) Then
          x = b(i1, i1)
          b(i1, i1) = b(i, i)
          b(i, i) = x
          For k = 1 To m
            x = w(k, i1)
            w(k, i1) = w(k, i)
            w(k, i) = x
          Next k
        End If
      Next i
'C
      For i = 1 To m
        c = 0
        For k = 1 To m
          c = c + w(k, i)
        Next k
        If (c < 0#) Then
          For k = 1 To m
             w(k, i) = -w(k, i)
          Next
        End If
      Next i
'      Return
      End Sub

Private Sub INV(a(), m, L, DET)
'C
'C     INVERSE AND DETERMINANT OF A MATRIX
'C      <INV>                            '98.4 '99.9
'C
'      IMPLICIT REAL*8(A-H,O-Z)
      'DIMENSION A(L,L),B(L,L)
       ReDim b(L, L)
      E = 0.00000001
      DET = 1
      For i = 1 To m
        For j = 1 To m
          b(i, j) = 0
          If (i = j) Then b(i, j) = 1
        Next j
      Next i
'C
      For k = 1 To m
        K1 = k
        c = Abs(a(k, k))
        For i = k + 1 To m
          If (Abs(a(i, k)) > c) Then
            K1 = i
            c = Abs(a(i, k))
          End If
         Next i
        If (c < E) Then
          MsgBox "MATRIX IS SINGULAR !!!! 逆行列を求めることができません"
         End
        End If
        If Not (k = K1) Then
          DET = -DET
          For j = 1 To m
            x = a(K1, j)
            a(K1, j) = a(k, j)
            a(k, j) = x
            x = b(K1, j)
            b(K1, j) = b(k, j)
            b(k, j) = x
          Next j
        End If
'C
        DET = DET * a(k, k)
        p = a(k, k)
        For j = 1 To m
          a(k, j) = a(k, j) / p
          b(k, j) = b(k, j) / p
        Next j
        For i = 1 To m
          If Not (i = k) Then
            q = a(i, k)
            For j = 1 To m
              a(i, j) = a(i, j) - q * a(k, j)
              b(i, j) = b(i, j) - q * b(k, j)
            Next j
          End If
        Next i
      Next
'C
      For i = 1 To m
        For j = 1 To m
          a(i, j) = b(i, j)
        Next j
      Next i
 '     Return
End Sub


Private Sub coloring(irow, jcolomn, icolor)
'34 水色 4 緑
        Cells(irow, jcolomn).Select
        With Selection.Interior
         .ColorIndex = icolor
         .Pattern = xlSolid
         .PatternColorIndex = xlAutomatic
         End With

End Sub

Public Sub 回転Oblimin()
Dim st$, j
   j = 1
   st$ = " Factor pattern matrix, Communality and Unique variance"
'st$ = " Orthomax rotated factor pattern matrix, Communality and Unique variance"
   Do
     j = j + 1
   Loop Until (Cells(j, 1) = st$) Or (j > 10000)
   If j > 10000 Then
     MsgBox (st$ + "の解を見つけることができませんでした。終了します。")
     End
   End If
   j = j + 1 '+ 1

'因子数を調べる
   k = 2
   Do
     k = k + 1
   Loop Until (Cells(j, k) = "" Or k > 25)
   If k > 25 Then
     MsgBox ("因子数がわかりませんでした。終了します。")
     End
   End If
   m_fact = Cells(j, k - 1)
'   MsgBox m_Fact

'変数の数を数える
   k = j
   Do
     k = k + 1
  Loop Until (Cells(k, 1) = "Variance" Or Cells(k, 1) = "" Or k - j > 200)
   If (k - j) > 200 Then
     MsgBox ("因子数がわかりませんでした。終了します。")
     End
   End If
  nv = k - j - 1

  ReDim a_fact(nv, m_fact), varnum(nv), varnames(nv) As String
   For jj = 1 To nv
     For ii = 1 To m_fact
       a_fact(jj, ii) = Cells(j + jj, ii + 2)
     Next ii
     varnum(jj) = Cells(jj + j, 1)
     varnames(jj) = Cells(jj + j, 2)
   Next jj
' 書き込み開始セルを見つける。、

   istart = 10000
   Do
     istart = istart - 1
   Loop Until (Cells(istart, 1) <> "" Or istart <= 1)
   If istart <= 1 Then
     MsgBox "書き込み開始位置を決めることができませんでした。"
     End
   End If
   
   istart = istart + 1
 '  For jj = 1 To nv
 '    Cells(jj + istart + 2, 1) = varnum(jj)
 '    Cells(jj + istart + 2, 2) = varnames(jj)
 '  Next jj
   Cells(istart, 1).Activate

   Call oblimin(nv, m_fact, a_fact, varnum, varnames)
   Call sorting(istart + 5, 3, m_fact + 3)

End Sub


Private Sub oblimin(nv, m_fact, a_fact(), varnum(), varnames() As String)
ReDim PHI(m_fact, m_fact)
Dim kaiser, ak1, ak2, ak3, ak4
Dim i, j, w
Dim it, iconv As Boolean
ReDim sp(nv, m_fact)
Dim iout(4)
  For i = 0 To 4
      iout(i) = 1
  Next i
    
ak1 = -0
ak2 = m_fact
ak3 = 0
ak4 = -m_fact
kaiser = 1

Call DGQ1(nv, m_fact, a_fact, sp, PHI, kaiser, ak1, ak2, ak3, ak4, it, iconv)

'-------------------------------------------------------------------------
ii = ActiveCell.Row

'Call Title(aa$, 1, 1)
'Cells(2, 1) = "Ogasawara,H.(2000). ROEF Version 1.0 User's guide, http://www.res.otaru-uc.ac.jp/~hogasa/rosef10.html"
'(Cells(3, 1) = "の DGQ1をexcel vba に移植した。 by khori. 2001/9/13"
'-------------------------------------------------------------------------
Dim m
m = m_fact
If Not (iconv) Then
  Cells(ii, 1) = "direct oblimin (δ=0)は収束しませんでした。"
  Cells(ii + 1, 1) = "反復数="
  Cells(ii + 1, 2) = it
  Cells(ii + 2, 1) = "収束基準="
  Cells(ii + 2, 2) = eps
  Exit Sub 'End
  ii = ii + 2
End If
'因子パターン行列の表示
'Cells(ii, 1) = "反復数="
'Cells(ii, 4) = it
'Cells(ii, 6) = "収束基準="
'Cells(ii, 8) = eps
ii = ii + 1
    Cells(ii + 1, 1) = "因子パタン行列 (direct oblimin δ=0 )"
    Call cellname(ii + 1, 1, "_direct_oblimin因子パタン")
    
    ii = ii + 3
 '   For i = 1 To nv
 '     Cells(ii + i, 1) = Right$("   " + Str(i), 3) + ": " + varname(i)
 '   Next i
    For j = 1 To m
      Cells(ii, j + 2) = j
    Next j
    
    For i = 1 To nv
        Cells(ii + i, 1) = varnum(i)
        Cells(ii + i, 2) = varnames(i)
      For j = 1 To m
        
        Cells(ii + i, j + 2) = a_fact(i, j)
        Cells(ii + i, j + 2).NumberFormatLocal = "###0.000_ "
        If (Abs(a_fact(i, j)) > 0.4) Then Call coloring(ii + i, j + 2, 34)
       Next j
    Next i
    Range(Cells(ii + nv, 1), Cells(ii + nv, m + 1 + 1)).Select
    With Selection.Borders(xlEdgeBottom)
        .LineStyle = xlDouble
        .Weight = xlThick
        .ColorIndex = xlAutomatic
    End With
    Cells(ii, 1).Activate
    ii = ii + nv + 1
    Cells(ii, 1) = "寄与"
    For j = 1 To m
    w = 0
    For i = 1 To nv
        w = w + a_fact(i, j) ^ 2
    Next i
      Cells(ii, j + 1 + 1) = w ', "#0.0000")
    Next j
    Range(Cells(ii, 2 + 1), Cells(ii, m + 1 + 1)).NumberFormatLocal = "###0.000_ "
   ii = ii + 2
'End If

'回転後因子間相関行列の表示
If iout(4) = 1 Then
    Cells(ii, 1) = "因子間相関 ( direct oblimin δ=0 )"
    Call cellname(CInt(ii), 1, "_direct_oblimin因子間相関")
    
    ii = ii + 1
    For j = 1 To m
      Cells(ii, j + 2) = j
    Next j
    
    For i = 1 To m
      Cells(ii + i, 1) = " "
      Cells(ii + i, 2) = i
      For j = 1 To m
        Cells(ii + i, j + 2) = PHI(i, j)
        'Cells(ii + i, j + 2).NumberFormatLocal = "###0.000_ "
        Call corrcolor(ii + i, j + 2)
      Next j
    Next i
End If
ii = ii + m + 2

'因子構造行列の表示
    Cells(ii, 1) = "因子構造行列(direct oblimin δ=0 )"
    Call cellname(CInt(ii), 1, "_direct_oblimin因子構造行列")
    
    ii = ii + 1
    For i = 1 To nv
      Cells(ii + i, 1) = i
    Next i
    For j = 1 To m
      Cells(ii, j + 2) = j
    Next j
    
    For i = 1 To nv
        Cells(ii + i, 1) = varnum(i)
        Cells(ii + i, 2) = varnames(i)
      For j = 1 To m
        Cells(ii + i, j + 1 + 1) = sp(i, j)
        Cells(ii + i, j + 2).NumberFormatLocal = "###0.000_ "
        If (Abs(sp(i, j)) >= 0.4) Then Call coloring(ii + i, j + 2, 34)
      Next j
    Next i
End Sub


Private Sub DGQ1(m, mf, D(), DR(), PH(), kaiser, ak1, ak2, ak3, ak4, it, iconv As Boolean)
                '変数数、因子数、因子パタン、因子構造、因子間相関
      Dim imax, i, j, am, i1, i2, j1
      Dim R2, RS, s2, R4, R3S, R2S2, RS3, x2, X2R2, X2RS
      Dim x, k, b1, c1, d1, e1, b2, c2, d2, e2, b3, c3, d3, e3, b4, c4, d4, e4
      Dim a, b, c, d0, nr, f1, f2, f3, u1, u2, t1, t2, iw1
    '  IMPLICIT REAL*8(A-H,O-Z)
     ' Dim D(IR, IC),DR(IR, IC), PH(IC, IC),
      ReDim DOLD(m, mf), CO(4), BN(3), CM(m)
      iconv = True
      eps = 0.000001
      eps = 0.00001
      imax = 200
      am = m
      For i = 1 To m
        CM(i) = 0
        For j = 1 To mf
          CM(i) = CM(i) + D(i, j) ^ 2
        Next j
        If (kaiser = 1) Then
        For j = 1 To mf
          D(i, j) = D(i, j) / Sqr(CM(i))
        Next j
        End If
      Next i
 '     WRITE(6,109)(CM(I),I= 1 to M)
 ' 109 FORMAT(/,'COMM=',12F6.3)
'C
      For i = 1 To mf
        For j = 1 To m
          DR(j, i) = 0
        Next j
        For j = 1 To mf
          PH(i, j) = 0
          If (i = j) Then PH(i, j) = 1
      Next j, i
      For i = 1 To m
        For j = 1 To mf
          DOLD(i, j) = 10000
      Next j, i
      it = 0
10    it = it + 1
      For i1 = 1 To mf
        For i2 = 1 To mf
        If (i1 = i2) Then GoTo 4
  'C
        R2 = 0
        RS = 0
        s2 = 0
        R4 = 0
        R3S = 0
        R2S2 = 0
        RS3 = 0
        x2 = 0
        X2R2 = 0
        X2RS = 0
  'C
          For j = 1 To m
            x = 0
            For k = 1 To mf
              x = x + D(j, k) ^ 2
            Next k
            R2 = R2 + D(j, i1) ^ 2
            RS = RS + D(j, i1) * D(j, i2)
            s2 = s2 + D(j, i2) ^ 2
            R4 = R4 + D(j, i1) ^ 4
            R3S = R3S + D(j, i1) ^ 3 * D(j, i2)
            R2S2 = R2S2 + D(j, i1) ^ 2 * D(j, i2) ^ 2
            RS3 = RS3 + D(j, i1) * D(j, i2) ^ 3
            x2 = x2 + x
            X2R2 = X2R2 + x * D(j, i1) ^ 2
            X2RS = X2RS + x * D(j, i1) * D(j, i2)
          Next j
  'C
        b1 = 4# * PH(i1, i2) * x2 * R2 - 4# * x2 * RS
        c1 = 4# * x2 * R2 + 4# * PH(i1, i2) ^ 2 * R2 ^ 2 - 8# * PH(i1, i2) * R2 * RS + 4# * RS ^ 2
        d1 = 8# * PH(i1, i2) * R2 ^ 2 - 8# * R2 * RS
        e1 = 4# * R2 ^ 2
  'C
        b2 = 4# * PH(i1, i2) * X2R2 - 4# * X2RS
        c2 = 4# * X2R2 + 4# * PH(i1, i2) ^ 2 * R4 - 8# * PH(i1, i2) * R3S + 4# * R2S2
        d2 = 8# * PH(i1, i2) * R4 - 8# * R3S
        e2 = 4# * R4
  'C
        b3 = 4# * PH(i1, i2) * R2 ^ 2 - 4# * RS * s2
        c3 = (2# + 4# * PH(i1, i2) ^ 2) * R2 ^ 2 + 2# * R2 * s2 + 4# * RS ^ 2
        d3 = 4# * PH(i1, i2) * R2 ^ 2 - 4# * RS * R2
        e3 = 2# * R2 ^ 2
  'C
        b4 = 4# * PH(i1, i2) * R4 - 4# * RS3
        c4 = (2# + 4# * PH(i1, i2) ^ 2) * R4 + 6# * R2S2
        d4 = 4# * PH(i1, i2) * R4 - 4# * R3S
        e4 = 2# * R4
  'C
        d0 = ak1 * b1 + ak2 * b2 + ak3 * b3 + ak4 * b4
        c = ak1 * c1 + ak2 * c2 + ak3 * c3 + ak4 * c4
        b = ak1 * d1 + ak2 * d2 + ak3 * d3 + ak4 * d4
        a = ak1 * e1 + ak2 * e2 + ak3 * e3 + ak4 * e4
  'C
        CO(1) = 4# * a
        CO(2) = 3# * b
        CO(3) = 2# * c
        CO(4) = d0
  'C
  'C      WRITE(6,500)CO
  'C  500 FORMAT(//'CO=',4D12.4)
  'C
        Call CUBNE(CO, BN, nr)
  'C
        If (nr = 0) Then
           MsgBox ("CANNOT SOLVE THE CUBIC EQUATION !!!!" _
           + Chr(10) + Str(CO(1)) _
           + Chr(10) + Str(CO(2)) _
           + Chr(10) + Str(CO(3)) _
           + Chr(10) + Str(CO(4)))
  '      WRITE(6,100)(CO(I),I= 1 to 4)
  '  100 FORMAT(//'CANNOT SOLVE THE CUBIC EQUATION !!!!'/4D15.7)
  '      STOP
         End
        End If
     ''   IF(NR = 1) then go to 6
        If (nr >= 2) Then
          f1 = a * BN(1) ^ 4 + b * BN(1) ^ 3 + c * BN(1) ^ 2 + d0 * BN(1)
          f2 = a * BN(2) ^ 4 + b * BN(2) ^ 3 + c * BN(2) ^ 2 + d0 * BN(2)
          If (f2 < f1) Then BN(1) = BN(2)
          If (f2 < f1) Then f1 = f2
          If (nr = 2) Then GoTo 6
          f3 = a * BN(3) ^ 4 + b * BN(3) ^ 3 + c * BN(3) ^ 2 + d0 * BN(3)
          If (f3 < f1) Then BN(1) = BN(3)
        End If
  'C
6       'Next
        u2 = BN(1)
        u1 = Sqr(1# + 2# * PH(i1, i2) * u2 + u2 ^ 2)
        t1 = 1# / u1
        t2 = u2 / u1
        For j = 1 To m
          D(j, i2) = -u2 * D(j, i1) + D(j, i2)
          D(j, i1) = u1 * D(j, i1)
        Next j
        For j = 1 To mf
          If (j <> i1) Then PH(i1, j) = PH(i1, j) * t1 + PH(i2, j) * t2
          If (j <> i1) Then PH(j, i1) = PH(i1, j)
        Next j
  'C
4
        Next i2
        Next i1
'C
      For i = 1 To m
        For j = 1 To mf
          If (Abs(D(i, j) - DOLD(i, j)) >= eps) Then
          If (it > imax) Then
            iconv = False
            eps = Abs(D(i, j) - DOLD(i, j)) '*************
'            Exit Sub
          End If
          For i1 = 1 To m
            For j1 = 1 To mf
              DOLD(i1, j1) = D(i1, j1)
            Next j1
          Next i1
      GoTo 10
      End If
     Next j, i
'C
'      WRITE(6,101)KAISER,AK1,AK2,AK3,AK4,IT
'  101 FORMAT(//'DIRECT OBLIQUE METHOD WITH GENERAL QUARTIC CRITERIA'
'     &//'KAISER=',I2,' : AK1,AK2,AK3,AK2=',4F5.1//
'     &  'NO. OF ITERATIONS=',I4)
      If (kaiser = 1) Then
        For i = 1 To m
         For j = 1 To mf
           D(i, j) = D(i, j) * Sqr(CM(i))
         Next j
        Next i
      End If
      For i = 1 To m
        For j = 1 To mf
          DR(i, j) = 0
          For k = 1 To mf
             DR(i, j) = DR(i, j) + D(i, k) * PH(k, j)
          Next k
        Next j
      Next i
'      WRITE(6,103)
'  103 FORMAT(/,'ROTATED LOADINGS',/)
'      for 13 I= 1 to MF
'   13 WRITE(6,104)(D(J,I),J= 1 to M)
'  104 FORMAT(5X,12F6.3)
'      WRITE(6,105)
'  105 FORMAT(/'STRUCTURES'/)
'      for 16 I= 1 to MF
'         WRITE(6,104)(DR(J,I),J= 1 to M)
'      Next i
'      WRITE(6,107)
'  107 FORMAT(/'FACTOR CORRELATIONS'/)
'      For i = 1 To MF
'        WRITE(6,104)(PH(I,J),J= 1 to MF)
'      Next i
End Sub
      
'C
'C     NEWTON METHOD FOR CUBIC EQUATION
Private Sub CUBNE(CO(), BN(), m)

     Dim eps, imax, in0, dd, a1, a2, a3, h, ai, fh, x1, x2, f1, f2, x, gra, hes
     Dim d1, d2
'      IMPLICIT REAL*8(A-H,O-Z)
   '   DIMENSION CO(4), BN(3)
      eps = 0.000000000001
      imax = 200
      in0 = 0
'C
      If (CO(1) = 0#) Then
'        WRITE(6,100)(CO(I),I= 1 to 4)
'  100 FORMAT(20X,' EQUATION IS NOT CUBIC ',30('!'),/4D20.13)
 
        If (CO(2) = 0# And CO(3) <> 0#) Then
          m = 1
          BN(1) = -CO(4) / CO(3)
          Exit Sub
        End If
        If (CO(2) = 0# And CO(3) = 0#) Then
          m = 0
          Exit Sub
        End If
'C
        dd = CO(3) ^ 2 - 4# * CO(2) * CO(4)
        If (dd < 0#) Then m = 0
        If (dd < 0#) Then Exit Sub
        m = 2
        BN(1) = (-CO(3) + Sqr(dd)) / (2# * CO(2))
        BN(2) = (-CO(3) - Sqr(dd)) / (2# * CO(2))
        Exit Sub
      End If
'C
'C     INITIAL VALUE FOR NEWTON METHOD
      a1 = CO(2) / CO(1)
      a2 = CO(3) / CO(1)
      a3 = CO(4) / CO(1)
      dd = a1 ^ 2 - 3# * a2
      h = -a1 / 3#
      ai = h
      fh = h ^ 3 + a1 * h ^ 2 + a2 * h + a3
'C
      If (dd < 0#) Then GoTo 1
      If (dd = 0#) Then ai = h - fh
      If (dd = 0#) Then GoTo 1
'C
      x1 = (-a1 - Sqr(dd)) / 3#
      x2 = (-a1 + Sqr(dd)) / 3#
      f1 = x1 ^ 3 + a1 * x1 ^ 2 + a2 * x1 + a3
      f2 = x2 ^ 3 + a1 * x2 ^ 2 + a2 * x2 + a3
'      WRITE(6,102)F1,F2
'  102 FORMAT('F1,F2=',2D21.14)
      If (f1 = 0#) Then x = x1
      If (f1 = 0#) Then GoTo 2
      If (f2 = 0#) Then x = x2
      If (f2 = 0#) Then GoTo 2
      If (f1 * f2 > 0#) Then
        If (f1 > 0#) Then ai = x1 - Sqr(-2# * f1 / (6# * x1 + 2# * a1))
        If (f1 < 0#) Then ai = x2 + Sqr(-2# * f2 / (6# * x2 + 2# * a1))
      End If
1     x = ai
'C
'C     ITERATION
      in0 = 0
      Do
3     in0 = in0 + 1
      gra = x ^ 3 + a1 * x ^ 2 + a2 * x + a3
      hes = 3# * x ^ 2 + 2# * a1 * x + a2
'C      WRITE(6,101)in0,GRA,HES,X
'  101 FORMAT('in0=',I4,1X,' F=',D21.14,' H=',D21.14,' X=',D21.14)
      If (in0 > 10 And (Abs(gra) < eps Or in0 > imax)) Then Exit Do
      x = x - gra / hes
      Loop
      'GoTo 3
'C
2    ' Next
      If (in0 > imax) Then
        m = 0
'        WRITE(6,101)in0,GRA,HES,X
        Exit Sub
      End If
      
      BN(1) = x
      d1 = a1 + x
      d2 = a2 + x * d1
      dd = d1 ^ 2 - 4# * d2
      If (dd < 0#) Then m = 1
      If (dd < 0#) Then Exit Sub
      m = 3
      BN(2) = (-d1 + Sqr(dd)) / 2#
      BN(3) = (-d1 - Sqr(dd)) / 2#
     
End Sub

Private Sub corrcolor(i, j)
     Cells(i, j).Select
     b = Abs(Cells(i, j))
     Selection.NumberFormatLocal = "##0.000_ "
     If b >= 0.7 And b < 1# Then
       Call coloring(i, j, 26)  'pink
       Else
        If b < 0.7 And b >= 0.5 Then
            Call coloring(i, j, 27)  'usumeno yellow
         Else
           If b < 0.5 And b >= 0.3 Then
                   Call coloring(i, j, 34)
            End If
        End If
     End If

End Sub

Private Sub cellname(i As Integer, j As Integer, x As String)
   Dim rc As Variant, rname As Variant
   rc = Cells(i, j).Address(ReferenceStyle:=xlR1C1) '"R2C5
'   rname = "=" + ActiveSheet.Name + "!" + rc
      rname = "=" + ActiveSheet.Name + "!" + rc
   If Len(ActiveSheet.Name) > 4 Then
     rc = ActiveSheet.Name + "_" + x
'     rc = x
   
   Else
     rc = ActiveSheet.Name + "__" + x
   End If
'   rc = "分析2__4xxxxxxxxxxxxxxxxxxxxxxxxx"
    ActiveWorkbook.Names.Add Name:=rc, RefersToR1C1:=rname
    'ActiveWorkbook.Names.Add Name:="myName", RefersToR1C1:= "=Sheet1!R1C1"

End Sub

Private Sub JACOBI(R(), n, ve(), eta)
       'Sub JACOBI(R(), n, e(), ve(), eta)
       '
       '固有値問題の解法（ヤコビ法）
'JACOBI:
      Dim i, j, k, AMX, L, NJ, NK
      Dim rr, rj, RK, rp, rm, w
      Dim c2, s2, cc, SC, VJ, VK
      ReDim E(n)  ' add
If n = 1 Then
      E(1) = R(1, 1)
      ve(1, 1) = 1
      Exit Sub
End If

      For i = 1 To n
        For j = 1 To n
            ve(i, j) = 0
        Next j
          ve(i, i) = 1
      Next i
       ''
    GoSub SEARCH_MAX      '非対角要素中での最大要素の探索
       '収束チェック
    While AMX > eta
         GoSub ROTATE       '最大要素をOにするように回転
         GoSub SEARCH_MAX
    Wend

       '対角要素が固有値
    For i = 1 To n
          E(i) = R(i, i)
    Next i

       '固有値の大きさの順に整列
    For i = 1 To n - 1
          L = i
           For j = i + 1 To n
           If E(L) < E(j) Then L = j
          Next j
          If i <> L Then
           Call SWAP(E(i), E(L))
           For k = 1 To n
            Call SWAP(ve(k, i), ve(k, L))
           Next k
          End If
    Next i
    For i = 1 To n
         R(i, i) = E(i)
    Next i
    Exit Sub

SEARCH_MAX:
    AMX = R(2, 1): NJ = 2: NK = 1
       '最大要素の探索
      For i = 2 To n
        For j = 1 To i - 1
           If Abs(R(i, j)) > AMX Then
            AMX = Abs(R(i, j))
            NJ = i: NK = j
           End If
        Next j
      Next i
      Return

ROTATE:
       ''R(NJ,NK) を０にする回転角 C の計算
    rr = R(NJ, NK)
    rj = R(NJ, NJ)
    RK = R(NK, NK)
    rp = (rj + RK) / 2
    rm = (rj - RK) / 2
    w = Sqr(rr * rr + rm * rm)
    c2 = Abs(rm) / w
    s2 = -Sgn(rm) * rr / w
       If rm = 0 Then s2 = -1
       'COS(c), SIN(c) の値
    cc = Sqr((1 + c2) / 2)
    SC = s2 / (2 * cc)
       '
       '回転後の値の計算
      For i = 1 To n
          rj = R(i, NJ)
          RK = R(i, NK)
          R(i, NJ) = rj * cc - RK * SC
          R(NJ, i) = R(i, NJ)
          R(i, NK) = rj * SC + RK * cc
          R(NK, i) = R(i, NK)
        
          VJ = ve(i, NJ)
          VK = ve(i, NK)
          ve(i, NJ) = VJ * cc - VK * SC
          ve(i, NK) = VJ * SC + VK * cc
      Next i
    w = rm * c2 - rr * s2
    R(NJ, NJ) = rp + w
    R(NK, NK) = rp - w
    R(NJ, NK) = 0
    R(NK, NJ) = 0
      Return
End Sub

Private Sub SWAP(a As Variant, b As Variant)
      Dim work

      work = a
      a = b
      b = work
End Sub

Private Sub sorting(istart, jstart, jend)
igeta = istart - 1
For i = 1 To nv
  amax = Abs(Cells(igeta + i, jstart))
  imax = jstart
  For j = jstart To jstart + m_fact - 1
    dum = Abs(Cells(igeta + i, j))
    If dum > amax Then
       amax = dum
       imax = j
    End If
  Next j
  Cells(igeta + i, jend + 2) = imax - jstart + 1
  Cells(igeta + i, jend + 3) = amax
Next i

Range(Cells(igeta + 1, 1), Cells(igeta + nv, jend + 3)).Select
Selection.Sort Key1:=Range(Cells(igeta + 1, jend + 2).Address), Order1:=xlAscending, _
       Key2:=Range(Cells(igeta + 1, jend + 3).Address), Order2:=xlDescending, _
       Header:=xlNo, OrderCustom:=1, MatchCase:=False _
       , Orientation:=xlTopToBottom, SortMethod:=xlPinYin
Cells(igeta, 1).Select

'作業領域の削除
Range(Cells(igeta + 1, jend + 2), Cells(igeta + nv, jend + 3)).Delete

End Sub