home *** CD-ROM | disk | FTP | other *** search
/ The Fred Fish Collection 1.5 / ffcollection-1-5-1992-11.iso / ff_progs / libs / matlab.lzh / MATLAB / MATLAB.LZH / Source / MatLab / MAGIC.FOR < prev    next >
Encoding:
Text File  |  1991-04-13  |  1.8 KB  |  82 lines
  1.       SUBROUTINE MAGIC (A, LDA, N)
  2.       IMPLICIT NONE
  3. C
  4. C ALGORITHMS FOR MAGIC SQUARES TAKEN FROM MATHEMATICAL RECREATIONS AND ESSAYS,
  5. C   12TH ED., BY W. W. ROUSE BALL AND H. S. M. COXETER
  6. C
  7.       INTEGER LDA, N
  8.       DOUBLE PRECISION A(LDA,N)
  9. C
  10.       INTEGER I, I1, IM, J, J1, JM, K, M, M1, M2, MM
  11.       DOUBLE PRECISION T
  12. C
  13. C
  14.       IF (MOD (N, 4).EQ.0) GO TO 100
  15.       IF (MOD (N, 2).EQ.0) M = N/2
  16.       IF (MOD (N, 2).NE.0) M = N
  17. C
  18. C ***      ODD ORDER OR UPPER CORNER OF EVEN ORDER
  19.       DO 20 J = 1, M
  20.         DO 10 I = 1, M
  21.           A(I,J) = 0
  22. 10      CONTINUE
  23. 20    CONTINUE
  24.       I = 1
  25.       J = (M+1)/2
  26.       MM = M*M
  27.       DO 40 K = 1, MM
  28.         A(I,J) = K
  29.         I1 = I-1
  30.         J1 = J+1
  31.         IF (I1.LT.1) I1 = M
  32.         IF (J1.GT.M) J1 = 1
  33.         IF (IDINT (A(I1,J1)).EQ.0) GO TO 30
  34.         I1 = I+1
  35.         J1 = J
  36. 30      CONTINUE
  37.         I = I1
  38.         J = J1
  39. 40    CONTINUE
  40.       IF (MOD (N, 2).NE.0) RETURN
  41. C
  42. C ***      REST OF EVEN ORDER
  43.       T = M*M
  44.       DO 60 I = 1, M
  45.         DO 50 J = 1, M
  46.           IM = I+M
  47.           JM = J+M
  48.           A(I,JM) = A(I,J)+2*T
  49.           A(IM,J) = A(I,J)+3*T
  50.           A(IM,JM) = A(I,J)+T
  51. 50      CONTINUE
  52. 60    CONTINUE
  53.       M1 = (M-1)/2
  54.       IF (M1.EQ.0) RETURN
  55.       DO 70 J = 1, M1
  56.         CALL RSWAP (M, A(1,J), 1, A(M+1,J), 1)
  57. 70    CONTINUE
  58.       M1 = (M+1)/2
  59.       M2 = M1+M
  60.       CALL RSWAP (1, A(M1,1), 1, A(M2,1), 1)
  61.       CALL RSWAP (1, A(M1,M1), 1, A(M2,M1), 1)
  62.       M1 = N+1-(M-3)/2
  63.       IF (M1.GT.N) RETURN
  64.       DO 80 J = M1, N
  65.         CALL RSWAP (M, A(1,J), 1, A(M+1,J), 1)
  66. 80    CONTINUE
  67.       RETURN
  68. C
  69. C ***      DOUBLE EVEN ORDER
  70. 100   CONTINUE
  71.       K = 1
  72.       DO 120 I = 1, N
  73.         DO 110 J = 1, N
  74.           A(I,J) = K
  75.           IF (MOD (I, 4)/2.EQ.MOD (J, 4)/2) A(I,J) = N*N+1-K
  76.           K = K+1
  77. 110     CONTINUE
  78. 120   CONTINUE
  79. C
  80.       RETURN
  81.       END
  82.