Pull to refresh

Comments 36

Прям сценарий для новой ветки HL1: молодой ученый писал программу для установки по изучению магнитного резонанса найденных обломков астероида, но ошибся в 2х строчках и в результате открылся портал в параллельный мир. Ну и далее мутанты, аномалии, пришельцы злые, тупые солдафоны и гречка по талонам!

Сделайте для потомков полезную вещь: напишите конвертер с Fortran 90 на C или C++.

Для Fortran 77 такое есть, а для Fortran 90 ещё нет.

Это возможно только в не имеющем большого практического значения грубом приближении, так как в C и C++, в отличие от Фортрана, не определён порядок выполнения арифметических операций в выражении.

Разве? Результат математических операций зависит от порядка их написания в этих языках. Скобки также есть.

В C и C++ порядок выполнения арифметических операций не определён, если алгебраически он не важен. Например, в выражении a+b+c компилятор Си вправе сначала сложить a и b, потом прибавить c, и также вправе сложить b и c, потом прибавить a. Компилятор Фортрана обязан складывать слева направо. Для плохо обусловленных численных задач (и при наличии побочных эффектов) это может быть важно.

Скобки в Си на это не влияют.

Operators can be regrouped according to the usual mathematical rules
only where the operators really are associative or commutative.

The additive operators + and - group left-to-right.

A parenthesized expression is a primary expression whose type and value
are identical to those of the enclosed expression. ... The parenthesized
expression can be used in exactly the same contexts as those where the
enclosed expression can be used, and with the same meaning, except as
otherwise indicated.

Ну да, вы ж сами процитировали самое основное:

Operators can be regrouped according to the usual mathematical rules only where the operators really are associative or commutative.

Поскольку арифметическое сложение – ассоциативная операция (в теории), то компилятор вправе её перегруппировывать.

Компилятор C/C++ вправе переупорядочивать операции любым образом, который в принципе сохраняет результат. Может даже вообще выкидывать – например, из a+b-a может сделать b.

А левая ассоциативность относится к логическому порядку понимания семантики выражения (в C++ для сложных объектов выполняется также и физически, но для простых численных типов может выполняться переупорядочивание).

Вот как написано в стандарте Си:

Order of evaluation of the operands of any C operator, including the order of evaluation of function arguments in a function‑call expression, and the order of evaluation of the subexpressions within any expression is unspecified...

А вот что написано про Фортран:

the expression is evaluated from left to right, according to the following precedence among operators...

Можно пример, когда использование скобок не решает проблему?

Да скобки вообще тут не причём и ничего не значат.

На практике большинство компиляторов Си в последнее время фактически использует фортрановский подход к последовательности арифметических операций, но стандарт языка этого не гарантирует.

А если говорить про собственно последовательность вычислений, а не только применения арифметических операций, то вот, например:

Hidden text
$ cat eval.c
#include <stdio.h>

int a () {
  puts ("a");
  return 1;
}

int b () {
  puts ("b");
  return 2;
}

int c () {
  puts ("c");
  return 3;
}

int main () {
  int s;
  s = a()*(b()*c());
  printf ("%d\n", s);
  return 0;
}

$ gcc eval.c
$ ./a.out
a
b
c
6

Это не к тому, что в Фортране кодогенерация устроена по-другому, а к тому, что скобки на неё не оказывают влияния.

Проверил на Java, результат такой же. Вы мне картину мира сломали :) Век живи, век учись! Как же тогда правильно сильно маленькие и сильно большие числа совместно обрабатывать правильно?

А если, как советуют, добавить компилятору флаги?

-fassociative-math, -ffast-math and -ffloat-store

Флагами, скорее всего, можно решить эту проблему, но это ж не выход в обсуждаемой ситуации. Чтобы понимать, какие флаги надо ставить в каждом конкретном случае, надо знать Фортран и Си на уровне, который сделает ненужным автоматический переводчик.

Есть какие-то рекомендации по написанию математических операций на сях? Скажем, имеем 2 малых величины m1, m2 и 2 больших величины, b1, b2.

Вместо (b1*b1)*(m1*m2) хотим (m1*b1)*(m2*b2) в предположении, что это уменьшит ошибку.

Обычно использование volatile отключает всю оптимизацию. Хотя это тоже не гарантировано стандартом.

Ну да, вы ж сами процитировали самое основное:

Operators can be regrouped according to the usual mathematical rules only where the operators really are associative or commutative.

Поскольку арифметическое сложение – ассоциативная операция (в теории), то компилятор вправе её перегруппировывать.

Интересно, что тут значит "really". Потому что сложение и умножение чисел с плавающей точкой как раз не ассоциативны. Так что я бы понял так, что перегруппировывать можно только операции над целыми числами.

К тому же для -fassociative-math в документации GCC явно сказано

Allow re-association of operands in series of floating-point operations. This violates the ISO C and C++ language standard by possibly changing computation result.

И как ни смешно, для Fortran это как раз может включиться по умолчанию:

For Fortran the option is automatically enabled when both -fno-signed-zeros and -fno-trapping-math are in effect.

До такой степени really операции над целыми числами тоже не ассоциативны, так как может происходить целое переполнение, которое, в свою очередь, может вызывать исключительную ситуацию или нет.

А как работают ключи у gcc/gfortran – это другой вопрос. С ключами по умолчанию входной язык не соответствует никакому стандарту.

Но выделенное явно подразумевает, что по стандартам C и C++ перегруппировывать операции с плавающей точкой нельзя.

И в MSVC:

By throwing the /fp:fast switch, you tell the compiler it should pretend that floats (and doubles) obey the rules of simple algebra (associativity and distributivity).

Do You Prefer Fast or Precise? - C++ Team Blog (microsoft.com)

То есть опять же без ключа перегруппировывать нельзя.

На всякий случай: я не говорю, что в примере из вашего комментария 28 мая в 19:40 (habr.com) определён порядок вычисления самих a(), b() и c(); он не определён, и я это хорошо знаю. Но только что скобки там причём, и гарантируют, что сначала будет вычислено `b() * c()`, а потом a() * результат.

Но выделенное явно подразумевает, что по стандартам C и C++ перегруппировывать операции с плавающей точкой нельзя.

... по мнению авторов GCC. В самом стандарте ничего такого не говорится. Стандарт Си вообще очень двусмысленно написан во многих моментах.

И в MSVC

Факт в том, что в последние годы компиляторы Си не переупорядочивают арифметические операции, так как от этого больше геморроя, чем пользы. Но стандарт ничего такого не требует. Поэтому и скобки ничего не гарантируют, хотя в обычных реализациях влияют.

До такой степени really операции над целыми числами тоже не ассоциативны, так как может происходить целое переполнение, которое, в свою очередь, может вызывать исключительную ситуацию или нет.

Я наконец нашёл исходную цитату, не знаю, почему в прошлый раз не получилось. И оказывается, там как раз говорится, что если целочисленное переполнение вызывает исключение, то перегруппировывать нельзя. [expr.pre] (eel.is)

[Note 4

The implementation can regroup operators according to the usual mathematical rules only where the operators really are associative or commutative.

For example, in the following fragment

int a, b;

/* ... */

a = a + 32760 + b + 5;

the expression statement behaves exactly the same as

a = (((a + 32760) + b) + 5);

due to the associativity and precedence of these operators.

Thus, the result of the sum (a + 32760) is next added to b, and that result is then added to 5 which results in the value assigned to a. 

On a machine in which overflows produce an exception and in which the range of values representable by an int is [-32768, +32767], the implementation cannot rewrite this expression as

a = ((a + b) + 32765);

...

However on a machine in which overflows do not produce an exception and in which the results of overflows are reversible, the above expression statement can be rewritten by the implementation in any of the above ways because the same result will occur.

 — end note]

Ключевое выделено. Так что чуть поправлюсь: перегруппировывать операции с плавающей точкой тоже можно, но только в тех случаях, когда это не может изменить результата (например, умножение на 0, когда можно доказать, что нигде нет бесконечностей или NaN).

Если вы найдёте утверждения разработчиков стандарта или хотя бы разработчиков компиляторов, что перегруппировывать с изменением результата по стандарту можно, готов изменить своё мнение.

Ну или можно сказать, что это Note и технически не обязательно к исполнению (non-normative) Element: Non-normative Note (iso.org)

Конвертер может вынести a+b в локальную переменную. Это разве не исправит проблему?

Гарантии нет. Не говоря даже о том, как это будет смотреться. Проще сразу на ассемблер.

Написание эффективного транслятора на С++ почти невозможно, ибо

  1. Fortran язык с исключительно гадкой грамматикой например program pz do do i=1,10 print *,"AAAAAAAAAAAA" enddo end program

  • валидно , а ваш лексический анализатор начинает сходить с ума.

  1. например, умножение матриц не ассоциативно с тз вычислительных затрат. Те. конвертер либо перегружает * и дает просадку по производительности, либо получаем трудно модифицируемый кусок кода будто через обсфукатор пропускали..

Фортран, тем самым, борется за выживание.

Но сейчас так написать нельзя.

огромный кусок с методом вращений заменён на единственную строку с вызовомnp.linalg.eigh

Вот интересно. Я сам никогда массивными расчетами не занимался, и с фортраном сталкивался только в университетском курсе программирования (достаточно бессмысленном, надо сказать). Но неоднократно слышал мнение, что фортран в научных расчетах не умирает, поскольку для него написаны мегатонны библиотек, от которых никто не может отказаться. И что, вы хотите сказать, что в библиотеке f90 не было более или менее стандартного способа получить собственные числа и векторы?

Да почему, LAPACK или MKL содержат всё нужное. Главное, не забывать RTFM.

А раньше была либо мода, либо правила хорошего тона - каждый под свои задачи кодит сам. В итоге имеем на кафедре десятки версий одного и того же с непредсказуемым набором багов. На одном из своих предметов я перестал требовать написание, например, метода Рунге-Кутты-Мерсона, устав ежегодно видеть ужастики и под копирку слизанные друг у друга программы. Вроде бы учить на готовых библиотеках, но с "плохими" тестовыми задачами получается легче и эффективнее по затратам времени.

Это да. У нас тут тоже есть некая программа для расчисления, так сказать, хода светил по небесной сфере. Телескоп наводить, куда надо. 250 килобайт вручную написанного некомментированного текста на f90 одним куском. У нее есть одна маленькая проблема: функции для работы с гражданским временем авторы почему-то решили написать сами. В те времена страна еще переходила каждый год с летнего времени на зимнее и обратно. И вот это они где-то там у себя захардкодили. Потом случился 2011 год... В общем, с тех пор это все обросло гроздьями внешних скриптов, которые с разной степенью успешности решают проблему. А в этом году вроде бы опять было предложение ввести переход на летнее время. Я как-то пытался с ней разобраться, но быстро сломался, заплакал и убежал.

Обычно для таких целей просто используют зафиксированное по отношению к UTC время на компьютере, как бы ни менялось гражданское.

Эта программа на Алголе написана до появления стандартных библиотек Фортрана.

В современном Фортране лучше не использовать специфические функции dabs, dsqrt и т.д. и константы вроде 1.0d0 – все стандартные функции полиморфны в отношении точности вещественных чисел. А так, если вы захотите в вашей программе поменять real*8 на real*10 или real*16, это не даст желаемого результата.

Это дурные привычки со студенчества, когда в наличии был жуткий Fortran Power Station 4.0. Постепенно переучиваюсь. Особенно удобно при работе с комплексными числами.

program NMR_rarefied
implicit none
real(8), dimension(:,:,:), allocatable :: dip, a1 ! dipolar interaction and additional
real(8), dimension(:,:), allocatable :: sx, sy, sz, sI, sp, sm, o, u, v, w, & ! spin matrices and additional
                                        x, y, z, r, & ! r = {x, y, z}, distances between magnetic ions (MI)
                                        J0, & ! exchange interactions
                                        ReH, ImH, b, q, q2, & ! Hamiltonian parts and additional
                                        d0, d1 ! for Kronecker product procedure
real(8), dimension(:), allocatable :: x1, y1, z1, & ! coordinates of MI
                                      RT2, RT3, RT4, RT1, sum_sqr, & ! for rotation method
                                      histp, hist ! partial and averaged histogram
logical, dimension(:), allocatable :: mask1
real(8), dimension(1:3) :: e1, e2, e3, d ! lattice basis
real(8), dimension(1) :: i11, k11
real(8) S, eps, & ! spin, precision
        alpha, beta, pd, & ! parameters of exchange (A*exp(-B*r)) and dipolar interaction (pd = (hbar*gamma)^2)
        res, & ! histogram resolution
        s1, s2, c1, c2, & ! rotation parameters
        m, e, c
integer N, N1, N2, N3, NN, Nres, & ! total number of MI and their concentration
        i, j, k, i1, j1, k1, t, m3, ms, l, r1, ires, & ! counters
        f, fn, & ! 2S+1, (2S+1)^N
        nh ! number of histogram bars
logical noconv
        
!----------------------------------------------------
! read and set the initial parameters
!----------------------------------------------------
open(1, file="init.txt", form="formatted", action="read")
read(1,*) eps, S, N, N1, N2, N3, nh, alpha, beta, pd, &
          e1, e2, e3, res, Nres
close(1)

f = 2*S + 1; fn = f**N
call RANDOM_SEED()

!----------------------------------------------------
! memory alloc
!----------------------------------------------------
allocate(sx(1:f,1:f)); allocate(sy(1:f,1:f)); allocate(sz(1:f,1:f))
allocate(sI(1:f,1:f)); allocate(sp(1:f,1:f)); allocate(sm(1:f,1:f))
allocate(o(1:f,1:f)); allocate(v(1:f,1:f)); allocate(w(1:f,1:f))
allocate(dip(1:6,1:N,1:N)); allocate(J0(1:N,1:N)); allocate(u(1:N,1:N))
allocate(x(1:N,1:N)); allocate(y(1:N,1:N)); allocate(z(1:N,1:N))
allocate(r(1:N,1:N))
allocate(x1(1:N)); allocate(y1(1:N)); allocate(z1(1:N))
allocate(ReH(1:fn,1:fn)); allocate(ImH(1:fn,1:fn))
allocate(b(1:fn,1:fn)); allocate(q(1:fn,1:fn)); allocate(q2(1:fn,1:fn))
allocate(mask1(1:fn))
allocate(a1(1:N,1:f,1:f))
allocate(RT1(1:fn)); allocate(RT2(1:fn))
allocate(RT3(1:fn)); allocate(RT4(1:fn))
allocate(sum_sqr(1:fn))
allocate(histp(1:nh)); allocate(hist(1:nh))

hist = 0.0d0

open(1, file="MI_coords.txt", form="formatted", action="write")
open(2,file="E.txt", form="formatted", action="write")

!----------------------------------------------------
! spin matrices initialization
!----------------------------------------------------
sx = 0.0d0; sy = 0.0d0; sz = 0.0d0; sp = 0.0d0; sm = 0.0d0; sI = 0.0d0
do i = 1, f, 1
  sI(i,i) = 1.0d0
  sz(i,i) = S - (i-1)
  if(i <= f-1) sp(i,i+1) = dsqrt(dabs(S*(S+1.0d0) - real(i*(i+1))))
enddo
sm = transpose(sp)
sx = 0.5d0*(sp + sm)
sy = 0.5d0*(sm - sp)

!----------------------------------------------------
ires = 1
REALIZATION: do while(ires <= Nres)
!----------------------------------------------------
! generation of MI distribution in lattice
!----------------------------------------------------
  x = 0.0d0; y = 0.0d0; z = 0.0d0; r = 0.0d0
  COORDS: do
  ! MI coordinates
    do i = 1, N, 1
      call RANDOM_NUMBER(d)
      d(1) = floor(d(1)*N1); d(2) = floor(d(2)*N2); d(3) = floor(d(3)*N3)
      x1(i) = d(1)*e1(1) + d(2)*e2(1) + d(3)*e3(1)
      y1(i) = d(1)*e1(2) + d(2)*e2(2) + d(3)*e3(2)
      z1(i) = d(1)*e1(3) + d(2)*e2(3) + d(3)*e3(3)
    enddo
  
  ! MI distances
    do i = 1, N, 1
      do j = 1, N, 1
        if(i /= j)then
          x(i,j) = x1(i) - x1(j)
          y(i,j) = y1(i) - y1(j)
          z(i,j) = z1(i) - z1(j)
          ! full interaction
          r(i,j) = dsqrt( x(i,j)**2 + y(i,j)**2 + z(i,j)**2 )
        endif
      enddo
      r(i,i) = 1.0d0
    enddo

    if(minval(r) > 0.5d0)exit
  enddo COORDS

  ! write MI coordinated
  write(1,'(A,I8)') "Realization ", ires
  do i = 1, N, 1
    write(1,'(I6, 3E20.6E3)') i, x1(i), y1(i), z1(i)
  enddo
  write(1,*)
  
!----------------------------------------------------
! pair interaction evaluation
!----------------------------------------------------
  ! isotropic part
  J0 = alpha*dexp(-beta*r)
  dip(1,:,:) = J0 + pd * (1.0d0 - 3.0d0*x**2 / r**2) / r**3
  dip(2,:,:) = J0 + pd * (1.0d0 - 3.0d0*y**2 / r**2) / r**3
  dip(3,:,:) = J0 + pd * (1.0d0 - 3.0d0*z**2 / r**2) / r**3
  ! anisotropic part
  dip(4,:,:) = -3.0d0*pd*x*y / r**5
  dip(5,:,:) = -3.0d0*pd*x*z / r**5
  dip(6,:,:) = -3.0d0*pd*y*z / r**5

!----------------------------------------------------
! real part of the Hamiltonian matrix
!----------------------------------------------------
  o = sz; call p78()
  ReH = b
  u = dip(1,:,:); v = sx; w = sx; call p56()
  u = dip(2,:,:); v = sy; w = -sy; call p56()
  u = dip(3,:,:); v = sz; w = sz; call p56()
  u = dip(5,:,:); v = sx; w = sz; call p56()
  q = ReH; ReH = 0.0d0

!----------------------------------------------------
! imaginary part of the Hamiltonian matrix
!----------------------------------------------------
  u = dip(4,:,:); v = sx; w = sy; call p56()
  u = dip(6,:,:); v = sy; w = sz; call p56()
  ImH = ReH; ReH = q

!----------------------------------------------------
! Diagonalization of the Hamiltonian matrix
!----------------------------------------------------
  t = 0; q = 0.0d0; q2 = 0.0d0
  do i = 1, fn, 1
    q(i,i) = 1.0d0
  enddo

  ! sum of squares by strings
  sum_sqr = 0.0d0
  do i = 1, fn, 1
    do j = 1, fn, 1
      if(i /= j) sum_sqr(i) = sum_sqr(i) + (ReH(i,j)**2 + ImH(i,j)**2)
    enddo
  enddo
  noconv = .false.
  
  ! rotation iterations
  ROTATION: do
    i11 = maxloc(sum_sqr); i1 = i11(1)
    mask1 = .true.; mask1(i1) = .false.
    k11 = maxloc(ReH(i1,:)**2 + ImH(i1,:)**2, mask=mask1); k1 = k11(1)
    m = maxval(ReH(i1,:)**2 + ImH(i1,:)**2, mask=mask1)
  
    if(m < eps)exit 
    if(t == 100000)then
      noconv = .true.
      exit
    endif

    mask1(k1) = .false.
    sum_sqr = sum_sqr - sum(ReH(:,i1)**2 + ImH(:,i1)**2, mask=mask1)
    
    i = i1; k = k1
  
    ! rotation parameters
    e = dabs(ReH(i,i) - ReH(k,k)) / dsqrt( (ReH(i,i) - ReH(k,k))**2 + 4.0d0*(ReH(i,k)**2 + ImH(i,k)**2) )
    s1 = dsqrt( 0.5d0*(1.0d0 - e) ); c1 = dsqrt( 0.5d0*(1.0d0 + e) )
    if(ReH(i,i) /= ReH(k,k)) s1 = sign(s1, ReH(i,i) - ReH(k,k))
    if(ReH(i,k) /= 0.0d0) s1 = sign(s1, ReH(i,k))
  
    s2 = ImH(i,k) / dsqrt(ReH(i,k)**2 + ImH(i,k)**2)
    if(ReH(i,k) /= 0.0d0) s2 = sign(s2, ReH(i,k))
  
    c2 = s1 * dabs(ReH(i,k)) / dsqrt(ReH(i,k)**2 + ImH(i,k)**2)
    s1 = s1*s2
  
    ! make the rotations
    RT1 = c1*q(:,i) + c2*q(:,k) + s1*q2(:,k)
    RT2 = -c2*q(:,i) + c1*q(:,k) + s1*q2(:,i)
    RT3 = -s1*q(:,k) + c1*q2(:,i) + c2*q2(:,k)
    RT4 = -s1*q(:,i) - c2*q2(:,i) + c1*q2(:,k)
    
    q(:,i) = RT1; q(:,k) = RT2; q2(:,i) = RT3; q2(:,k) = RT4

    RT1 = c1*ReH(:,i) + c2*ReH(:,k) + s1*ImH(:,k)
    RT2 = -c2*ReH(:,i) + c1*ReH(:,k) + s1*ImH(:,i)
    RT3 = -s1*ReH(:,k) + c1*ImH(:,i) + c2*ImH(:,k)
    RT4 = -s1*ReH(:,i) - c2*ImH(:,i) + c1*ImH(:,k)

    ReH(:,i) = RT1; ReH(:,k) = RT2; ImH(:,i) = RT3; ImH(:,k) = RT4

    RT1 = c1*ReH(i,:) + c2*ReH(k,:) - s1*ImH(k,:)
    RT2 = -c2*ReH(i,:) + c1*ReH(k,:) - s1*ImH(i,:)
    RT3 = c1*ImH(i,:) + c2*ImH(k,:) + s1*ReH(k,:)
    RT4 = -c2*ImH(i,:) + c1*ImH(k,:) + s1*ReH(i,:)

    ReH(i,:) = RT1; ReH(k,:) = RT2; ImH(i,:) = RT3; ImH(k,:) = RT4
    
    t = t+1
  
    mask1 = .true.; mask1(i1) = .false.; mask1(k1) = .false.
    sum_sqr = sum_sqr + sum(ReH(:,i1)**2 + ImH(:,i1)**2, mask = mask1) &
              + sum(ReH(:,k1)**2 + ImH(:,k1)**2, mask = mask1)
   
    mask1 = .true.; mask1(i1) = .false.
    sum_sqr(i1) = sum(ReH(i1,:)**2 + ImH(i1,:)**2, mask = mask1)
  
    mask1 = .true.; mask1(k1) = .false.
    sum_sqr(k1) = sum(ReH(k1,:)**2 + ImH(k1,:)**2, mask = mask1)
  enddo ROTATION

!----------------------------------------------------
  if(noconv)then
    print *, "Realization does not converge"
  else
    print *,ires, ", Rotation finished at ", t, " iterations"
!----------------------------------------------------
! eigenvalues
!----------------------------------------------------
    do i = 1, fn, 1
      RT1(i) = ReH(i,i)
    enddo

    ! write energy
    write(2,'(A,I8)') "Realization ", ires
    do i = 1, fn, 1
      write(2,'(I6,E20.6E3)') i, RT1(i)
    enddo
    write(2,*)

!----------------------------------------------------
! evaluate probabilities and histogram amplitudes
!----------------------------------------------------
    histp = 0.0d0
    o = sx

    call p78()

    ReH = matmul(b, q); ImH = matmul(b, q2)
    d1 = transpose(q); b = transpose(q2)

    q = matmul(d1, ReH) + matmul(b, ImH)
    q2 = matmul(d1, ImH) - matmul(b, ReH)

    do i = 1, fn, 1
      do j = 1, fn, 1
        c = q(i,j)**2 + q2(i,j)**2
        k = floor(1.0d0 + dabs(RT1(i) - RT1(j)) / res)
        if(k <= nh) histp(k) = histp(k) + c
      enddo
    enddo

    hist = hist + histp
    ires = ires + 1  
  endif
enddo REALIZATION
!----------------------------------------------------

open(3,file="hist.txt", form="formatted", action="write")
do i = 1, nh, 1
  write(3,'(2E20.6E3)') i*res, hist(i) / Nres
enddo
close(1); close(2); close(3)

print *,"Histogram is written at hist.txt"

!----------------------------------------------------
! memory clear
!----------------------------------------------------
deallocate(sx); deallocate(sy); deallocate(sz); deallocate(sI)
deallocate(o); deallocate(v); deallocate(w)
deallocate(ReH); deallocate(ImH); deallocate(b); deallocate(q); deallocate(q2)
deallocate(mask1); deallocate(a1)
deallocate(RT1); deallocate(RT2); deallocate(RT3); deallocate(RT4); deallocate(sum_sqr)
deallocate(histp); deallocate(hist)
deallocate(dip); deallocate(J0); deallocate(sp); deallocate(sm)
deallocate(x1); deallocate(y1); deallocate(z1)
deallocate(x); deallocate(y); deallocate(z); deallocate(r); deallocate(d1)

contains
!----------------------------------------------------
! Kronecker product of one spin and N-1 identity
!----------------------------------------------------
subroutine p78()
  b = 0.0d0

  do k = 1, N, 1
    a1(k,:,:) = sI
  enddo
  
  do k = 1, N, 1
    a1(k,:,:) = o
    call p34()
    a1(k,:,:) = sI
    b = b + d1
  enddo
end subroutine p78

!----------------------------------------------------
! Kronecker product of two spin and N-2 identity
!----------------------------------------------------
subroutine p56()
  do k = 1, N, 1
    do t = 1, N, 1
      if(k /= t)then
        a1(k,:,:) = v; a1(t,:,:) = w
      endif
      
      call p34()
            
      a1(k,:,:) = sI; a1(t,:,:) = sI
      ReH = ReH + u(k,t)*d1
    enddo
  enddo
end subroutine p56

!----------------------------------------------------
! Kronecker product of N matrices
!----------------------------------------------------
subroutine p34()
  if(allocated(d0)) deallocate(d0)
  if(allocated(d1)) deallocate(d1)
  
  allocate(d0(1:f,1:f)); d0 = 0.0d0
  allocate(d1(1:f**2,1:f**2)); d1 = 0.0d0
  
  d0 = a1(1,:,:)
  
  do m3 = 2, N, 1
    NN = f**(m3-1)
    do i = 1, NN, 1
      do j = 1, NN, 1
        do ms = 1, f, 1
          do l = 1, f, 1
            r1 = (i-1)*f + ms; i1 = (j-1)*f + l
            d1(r1,i1) = d0(i,j)*a1(m3,ms,l)
          enddo
        enddo
      enddo
    enddo
    
    if(m3 < N)then
      deallocate(d0); allocate(d0(1:f**m3,1:f**m3)); d0 = 0.0d0
      d0 = d1
      deallocate(d1); allocate(d1(1:f**(m3+1),1:f**(m3+1))); d1 = 0.0d0
    endif
    
  enddo
  deallocate(d0)
end subroutine p34

end program NMR_rarefied

вы бы хоть расскрасили. Например это "CSS",

Примерно так может выглядеть сигнал в эксперименте. Только колебания
происходят намного быстрее - десятки и сотни МГц и даже выше.

Намного быстрее чего? У вас на графике нет масштаба ни по времени, ни по напряжению.

Написано же, типичное время спада - десятки микросекунд. А частота колебаний - десятки и сотни МГц. Периодов несущей на соответствующей картинке будет гораздо больше.

Истинно так! :)
Поэтому, например, в статьях у меня графики выглядят будто со сплошной чёрной заливкой. Или с мелкой рябью.

Рандомная картинка из кучи промоделированных

Простите, не заглядывал в комментарии пару дней :)

Sign up to leave a comment.

Articles