Интеллектуальное суммирование в FOR(Mula) TRAN(slation)

Что ж, сегодня я здесь с этим сомнением...

Я хочу написать это уравнение на Фортране:

введите здесь описание изображения

Конечно, я могу использовать «классический» подход и написать так:

 do i=1,N
        ac=0.0
        do j=i+1,M
            ac=ac+A(i,j)*B(j)
        enddo
        B(i)=(B(i)-ac)/A(i,i)
  enddo

Но так как я пишу на Фортране, то хотелось бы написать выражение, похожее на "больше похожее на оригинал" и, таким образом, компактное. Я думал примерно так:

forall(i=1:N,j=i+1:M)
    B(i)=(B(i)- <MySummationExpression(A(i,j)*B(j))> )/A(i,i)
endforall

Выражение, которое гораздо больше похоже на оригинал. Но правда в том, что мне трудно найти способ написать выражение суммирования простым и компактным способом. Конечно, я могу написать функцию "real function summation(<expression>,<lower bound>, <upper bound>)", но поскольку мы говорим о Фортране, я подумал, что должен быть простой (может быть, встроенный (?)) способ ее написания.

Итак, есть ли компактный способ написать это выражение, или я должен выбрать более уродливый способ (два явных цикла)?

EDIT: В реальном коде x представляет собой двумерный массив с одним решением в каждом столбце. Таким образом, использование встроенной функции sum, которая до сих пор кажется отличной идеей (как показано в ответе @alexander-vogt), приводит к почти такой же «компактности» кода:

do j=1,size(B,2)
        do i=nA,1,-1
            B(i,j)=(B(i,j)-sum(A(i,i+1:nA)*B(i+1:nA,j)))/A(i,i)
        enddo
 enddo

algorithm fortran fortran90 fortran95
person Alfonso Santiago    schedule 09.10.2013    source источник
comment
Надлежащее усилие, но более умный написанный код часто менее удобочитаем и может содержать бесплатные ошибки. Пишите код так, как вам удобнее.   -  person milancurcic    schedule 10.10.2013

Ответы (1)


arrow_upward
2
arrow_downward

Как насчет:

do i=1,N
  B(i) = (B(i) - sum( A(i+1:M,i) * B(i+1:M) )) / A(i,i)
enddo

(Обратите внимание, что я изменил индексацию матрицы A, Fortran является основным столбцом!)

Поскольку мы повторно используем B для результата, этот цикл должен выполняться в порядке возрастания. Я не уверен, что это возможно сделать с помощью forall (разве компилятор не должен выбирать, как действовать? - см. Fortran для всех ограничений).

Запись результата в новый вектор C не перезаписывает B и может выполняться в любом порядке:

forall ( i=1:N )
  C(i) = (B(i) - sum( A(i+1:M,i) * B(i+1:M) )) / A(i,i)
endforall
person Alexander Vogt    schedule 09.10.2013
comment
Я вижу, что вы изменили индексацию, но таким образом результат неверен. Я буду исследовать дальше, почему. - person Alfonso Santiago; 10.10.2013
comment
Что ж, просто меняя на A(i,i+1:M), я получаю те же результаты, что и с кодом в исходном вопросе... Матрица A должна быть транспонирована для работы с предоставленным мной решением ;-) - person Alexander Vogt; 10.10.2013