물리학자를 위한 포트란 90

Transcription

1 물리학자를 위한 FORTRAN90 이인호 한국표준과학연구원 2003년 2월 23일 출발 FORTRAN (FORmula TRANslation, 포트란)이란 무엇인가? 포트란은 엔지니어, 이공학자 및 기타 다른 과학적 알고리즘의 제작자나 사용자들을 위해 설계된 3세대 프로 그래밍 언어이다. 1950년대 말에 IBM의 John Backus에 의해 기획/설계되었으며, 매우 간결하고 엄격한 구문 형식을 가지고 있으며 배우기 쉬운 것이 특징이다. 포트란의 가장 유명한 두 개의 버전은 FORTRAN IV와 FORTRAN77을 꼽을 수 있다. FORTRAN IV는 1966년에 USASI 표준으로 승인되었으며, FORTRAN77은 1978년에 ANSI에 의해 승인된 버전이다. ISO와 ANSI에 의해 포트란의 표준으로 새로 인정받은 것을 FORTRAN90이라고 부르며, 1990년 초에 개발되 었다. 포트란90은 상당한 수준의 기능을 갖춘 현대화된 언어이다. 이것은 기본적으로 프로그래밍 언어 C++ 에서 지원하는 객체지향 (object-oriented) 프로그래밍 기능들과 병렬 배열 문법을 포트란77 언어에 추가한 것 이다. 오늘날에는, C (C++) 또는 자바 언어가 포트란을 대부분 대체하였지만, 그러나 아직도 (?) 포트란 사용자들 이 많이 있다. 그 이유는 전산 물리학 계산에서 사용해오는 전산 소프트웨어들이 다세대에 걸쳐서 계속 사용 되기 때문이다. 또한 누구나 쉽게 배울 수 있다는 언어의 특징이 있다. 특정 연구그룹에서 이미 검정을 거친 루틴들이 있다. 물론, 공공적인 공간에서도 이러한 프로그램들은 존재한다. 이들을 사용하기 위해서는 포트 란 언어를 알아두는 것이 필요하다는 것을 깨닫게 되는 것이다. 실제로 응용이 많이 되는 루틴들은 저 마다 오랜 역사를 가지고 있으며 포트란으로 구현되어 있는 것들이 많이 있다. 이러한 루틴들은 적극적으로 사용 되어야 우리가 원하는 더 많은 일들을 할 수 있을 것이다. 포트란을 사용하지 않는 사람들도 이미 개발된 포 트란 프로그램들을 읽을 수 있으면 자신들의 일에 어떠한 형태로 이든지 도움이 될 수 있다. FORTRAN90의 특징 1. 문법이 매우 간단하여 배우기가 쉽다. (고수준 언어이다; 본질적으로 기계보다 사람에게 더 가깝 게 설계되어 있다.) 2. 과학 기술 계산용 전문 언어이다. (미리 정의된 함수들 다수 보유, 113, 다양한 수치 계산 프로그래 밍 모듈 지원) 3. 수식과 여러 종류의 계산등을 간단히 공식화할 수 있다. (감각적인 간결성, 묵시적 상황 판단, 추상 적 자료형 지원) 4. 대상중심적이다. (동적인 메모리 관리, 모듈화의 장점을 살릴 수 있다, 확실한 데이터 처리와 보완 성 확보) 5. 기종간의 호환성이 있다. (이미 검증된 많은 포트란77 프로그램과 쉽게 연동된다.) (1 of 114) 오전 10:21:34

2 객체지향(대상중심)이란?: 특정한 응용을 위한 객체 (object)를 기술하는 높은 수준의 데이터 형, 데이 터 형 처리 방법들을 정의 함으로써 프로그램을 단독으로 작성하는 것 (개발 후에 있게 되는 수정/보 완 시에도 단독으로...). 프로그램 전체에서 다른 부분 사이가 독립되어 있기 때문에 이해하기 쉽고 확 장과 변경이 용이한 코드를 얻을 수 있는 장점이 있다. 객체는 데이타와 함수의 모음으로 표현된다. 아 래에 나타낸 것이 포트란90과 관련된 객체지향성이다. 이들은 모듈의 작성과 깊이 연관되어 있다. 객 체지향이라는 말보다는 대상중심이라는 말이 더 적합하다고 생각합니다.특정 대상을 전문으로 취급 하는 데이터와 부 프로그램들을 설계하고 필요에 따라서 개발 기간동안 쉽게 계속 추가하는 것을 목 표로 한다. 1. data abstraction -- user-defined types 2. data hiding -- private and public attributes 3. encapsulation -- modules and data hiding facilities 4. inheritance and extensibility -- generic procedures 5. polymorphism -- generic overloading 6. reusability -- modules 참고 웹자료: Introduction to object-oriented concepts using fortran90 현대적 의미로서의 포트란은 새로 만들어지는 코드는 포트란90으로 작성하며, 기존의 포트란77로 작성된 코 드들은 철저히 재 활용 (recycling)한다는 의미이다. 즉, 검정이 완전히 끝난 코드들 ( 다시 포트란90으로 재 작성할 필요는 없으며 그러한 일들은 많은 경우 무의미할 뿐만 아니라 오히려 위험할 수 있 다. 프로그램의 테스트와 안정성을 유지하는 것이 그만큼 어려운 것이라는 것을 의미한다. 현대의 컴퓨터 언 어를 기술할 때 다음의 사항들이 중점적으로 거론된다. reusability, maintainability, readability, reliability, testability. 본 웹 페이지에서는 포트란90의 강력한 기능인 포트란 모듈의 작성에 관해서 서술하기로 한다. 현존하는 포 트란77로 작성된 코드도 사용자의 편의를 위해서 또는 사용자의 목적에 알맞게 포트란90의 모듈로 재설계/ 작성하면 보다 일반적이고도 안전한 모듈이 될 수 있다. 이러한 일들은 모듈 개념을 적용하면 아주 쉽게 할 수 있다. 최소한, 현존하는 포트란77 코드를 읽을 수 있어야 한다. 소위 엔지니어링과 사이언스를 하면서 프 로그램의 신주류를 찾고 깊이 이해하여 활용하기란 매우 어렵고 또한 그러한 것들이 최고의 상황이 아닐지 도 모른다. 이러한 의미에서 검정된 언어의 우수한 기능들을 점진적으로 연습하는 것이 상당히 합리적이라 할 수 있다. John Backus; 포트란의 아버지 (2 of 114) 오전 10:21:34

3 왜 모듈인가? FORTRAN77에서는 찾아 볼 수 없는 기능이다. 대상중심 (객체지향 프로그래밍; object-oriented programming; "대상중심이라는 말이 더 어울림.) 프로그램에서만 나오는 개념이다. 모듈을 사용하면 변수와 변수 관련 함 수(서브루틴 포함의 일괄적이고 보다 안전한 관리를 조장한다. 사실 FORTRAN77은 소위 "Structured Program"을 지향했다. 쉽게 비유하면, 단 한 번에 거대한 집을 설계하는 것이다, 그 다음, 그 설계도에 맞추어 서 충실히 집의 구석구석을 열심히 만들어 나가는 것이다. 물론, 이러한 아이디어 자체가 나쁘다는 것이 결 코 아니다. 아주 이성적인 선택이며 분명히 좋은 방법임에 틀림없다. 하지만, 객체지향적 프로그래밍이 훨씬 프로그래밍 효율면에서 우위를 점한다는 것이 경험적으로 알려지기 시작했다. 소위 프로그램 단위의 '재사용 (reusability)'과 변형 (modification)이 아주 쉽게, 안전하게 진행될 수 있다. 프로 그램은 새로운 환경과 아이디어에 의해서 쉽게 변형되어져야 한다. 왜냐하면 프로그램은 우리의 아이디어 를 적용하고 테스트하는 방법이기 때문이다. 프로그래밍이 우리의 아이디어 창출에 장애가 되면 결코 안된 다. 이러한 의미에서 대상중심의 프로그래밍이 구조적 프로그래밍 보다 매우 유리하다. 특별한 아이디어를 빨리 구현하는 것이 얼마나 중요한가는 지난 컴퓨터 프로그램의 역사를 통해서 증명되었다. 속도는 차후에 개선할 수 있지만, 독창적인 아이디어의 발표는 그렇지 못하다. "Structured Program"에서는 작성하는 프로그램이 커질수록, 다양한 전체 프로그램의 시도가 있을수록, 여러 사람과 공동으로 일할수록, 시간이 지나서 버전이 올라 갈수록, 무슨 작업들이 어디에서 일어나는지, 어떠한 계산이 진행되는지 또한 변수들에 어떠한 변화들이 일어나는지 추적하기가 점점 어려워진다. 다시 말해서, 점점 더 암호화되고 복잡해지고 그에 따른 관리가 힘들어진다는 것이다. 따라서, 새로운 프로젝트를 위한 생 산성이 상당한 수준에서 위협받을 수 있다. 모듈을 사용하면 프로그램에 전체적인 변화가 분명히 있음에도 불구하고 구체적인 소스코드 상에서는 국소 적인 코드 변화로 종결되는 것을 조장한다. (위에서 기술한 내용은 다소 추상적이라고 평가받을 수 있다. 보 다 냉정하고 엄밀히 이야기하면 프로그래머의 습관에 관한 사항이다. 특히, 모듈을 사용하지 않는 사람들에 게는 대부분 이해가 되지 않을 뿐만 아니라 언급된 대부분의 항목에 대해서 하나 하나 변론할 수도 있다는 것 을 잘 알고 있다. 이것은 명백한 현실이다. 지난 세월, 사실, 세상의 프로그램들 모듈 없이도 잘 되어 왔다. 결 국, 위에서 언급한 모듈의 진정한 의미는 실제 모듈을 사용해 본 다음에 평가가 내려져야 마땅할 것이다. 모 듈들을 잘 설계하고 잘 사용하고 나면 그 유명한 "순서도" 조차도 특별한 의미와 필요가 없어진다면 믿어지 겠는가? 주 함수에는 정말 '말하듯' 주요 명령들만이 줄지어 나타나게 된다는 것이다. 이 때 부터는 정말 프로 그래머에게서 아이디어가 살아나도록 프로그램이 유도하는 수준이라고 하면 너무 심한 비약일까?) FORTRAN90 모듈은 무엇인가? 프로그래밍 언어 C++의 class와 비슷한 것이다. 포트란77의 common, include, block data를 대치할 수 있다. 물 론, 이것이 다는 아니다. 각종 형태의 데이터와 이들 데이터와 관련된 함수들의 연합체가 모듈이다. 다시 말 해서 객체를 나타내는 수단으로 모듈이 있을 수 있다. 여기서 데이터는 미리 정해진 사이즈의 배열일 필요가 없다. 배열의 크기를 정하는 함수를 모듈 속에 하나 만들어 두고 사용할 수 있다. FORTRAN77에서는 공통으로 사용하는 변수들을 모아 두고 사용하는 common 문이라는 것이 있었다. common 문은 FORTRAN77에서의 인수를 제외할 경우 유일한 global 변수의 접근 통로로 사용되었다. 소위 named common 문을 사용하여 몇 가지로 그 global 변수들을 편의상 분류하기도 했다. 하지만, 불안정한 프로 (3 of 114) 오전 10:21:34

4 그래밍의 이유가 된다고 하여 그것의 사용이 과도하게 느낄 정도로 엄청나게 천대와 탄압을 받아왔다. 하지 만, FORTRAN90에서는 보다 안전하게 모듈을 사용하여 common 문 그 고유의 장점을 살릴 수 있다. 또한 제 도적으로 보다 안전하게 (save, use, public, only, => 키워드를 사용함으로써) 프로그램할 수 있게 되었다. 굳이 변명하자면 common 문 속의 변수들을 사용하는 함수들은 거의 정해져있다. 따라서, 이들 변수들과 이 를 이용하는 함수들은 붙어서 다니는 것이 좋겠다. 이것이 모듈이다. 저장된 변수들의 접근이 명시될 수 가 있기 때문에 공통으로 사용되는 변수들의 사용에서 그 투명성이 증대된다. 여러 명이 프로그램을 개발하여 합칠 때에도 각 모듈들 간의 있을지도 모를 충돌을 제도적으로 방지할 수 있게 했다. 사실상 신경을 곤두세워 서 보아야 할 변수들을 키워드 public을 이용해 목록해 두도록 되어있다. 이어지는 프로그램 수정과 확장에서 도 그 탁월한 유용성이 보장된다. 자주 일어나는 일이지만, 특정 데이터를 생성하는 곳과 그 데이터를 참조하여 사용하는 곳이 멀리 (인수로 통해서 접근한다고 생각할 때) 떨어져 있으면, common 문을 사용하면 쉽게 필요로 하는 정보를 공유할 수 있 다. 구차하게 필요한 정보들을 전달하고 또 전달하는 귀찮은 작업들을 할 필요가 없다. 문제는 너무 많은 작 업들을 이런 식으로 하는 것은 결국 그러한 정보 관리상의 어려움을 유발하고 만다. 동적으로 크기가 변화하는 변수들도 포트란90에서 지원되기 때문에 정확히 그 크기와 순서가 일괄적으로 정 해져있는 포트란77에서의 common문 한계를 쉽게, 확실하게 극복할 수 있다. 사실, 이러한 위험성 때문에, 포 트란77에서는 include문을 사용하여 완벽하게 똑같은 common문들을 사용했었다. (예, common문을 필요로 하는 서브루틴마다, include 'file_common_statements'를 선언하고, 해당파일만 다루는 문제마다 바꿈으로서, 일괄적으로 관리할 수 있었다.) 어떻게 해야 모듈을 잘 만들고 잘 사용할 수 있는가? 모듈을 작성할 때는 save, public (또는 private)기능들을 확실히 사용하고, 모듈을 이용할 때는 이들을 꼭 확인 한다. 각각의 루틴을 작성할 때 implicit none을 항상 사용한다. 남이 만든 모듈과 내가 만든 프로그램 상에서 우연히 서로 다른 변수를 의미함에도 불구하고 변수 이름이 동등할 수 있다. 이 경우 그 숫자가 너무 많지 않 고 오해의 소지가 없다면 간단히 변수이름을 재명명함으로써 슬기롭게 문제를 피해갈수 있다. 이것 또한 포 트란90의 장점이다. 위에서 지적한 implicit none은 모든 변수들의 속성이 반드시 선언되어야 한다는 의미이 다. 즉, 뒤 따르는 줄에서 실수형, 문자형, 배열형, 정수형, 논리형 등이 구체적으로 정의되어야만 한다는 의미 이다. 그렇지 않은 변수들의 리스트를 컴파일 시 보여주고 더 이상의 컴파일을 진행하지 않는다. 포트란에서 는 컴파일할 때 보통 -C 옵션을 이용하면 배열의 한도를 넘는 배열 참조에 관련된 에러발생 위치를 실행할 때 보여주고 실행이 중단된다. 프로그램 개발 단계에서 필요한 옵션이다. USE stack, local_use_pop => pop pop가 특정 함수에서 국소적으로 사용될 때, stack 모듈 속의 pop는 local_pop로 재명명되었다. USE stack, ONLY: pos, local_pop=> pop 모듈 stack속의 pos와 pop만이 사용 가능하다. 다만, 이 경우 모듈 속 pop의 경우 local_pop로 이름이 다시 만들 어졌다. 모든 모듈들은 그것을 사용하는 프로그램 단위 보다 먼저 컴파일되어야 한다. 이렇게 되기 위해서는 makefile (4 of 114) 오전 10:21:34

5 에서 아래와 같은 프로그램들간의 의존성이 필요하기도 하다. 즉, 프로그램 ximage_lbfgs.f90은 모듈 timpose. f90를 사용하게 되어있다. 이 때 이 모듈은 다른 파일 속에 있는 것이며 이전에 만들어진 것이다. USE aaa USE bbb 모듈 bbb에서 public한 모든 변수, 배열, 또는 부 프로그램(function, subroutine)들은 모두 사용할 수 있다. USE ccc USE ddd, ONLY: qqq,ppp 모듈 ddd에서 사용하고자하는 변수, 배열, 또는 부 프로그램 (subroutine, function) 이름들을 적어둔다. 특정 모듈을 사용할 때, "USE 모듈이름은 아주 앞 쪽에 선언된다. 즉, 그 모듈을 사용하는 서브루틴에서 모든 변수들 보다 앞서서 선언된다. 결국, implicit none 보다 앞서서 그 바로 윗줄에서 정의된다. 즉, 컴파일할 때 먼 저 불러져서 컴파일이 될 수밖에 없다. USE lll, ONLY: sss 등과 같이 ONLY를 많이 사용할수록 프로그램의 readability는 당연히 높아지게 된다. 보다 더 확실하게 변수들을 사용한다는 것이다. 일반적으로 ONLY를 사 용하는 것이 사용하지 않는 것보다 유리하다. 구조적으로 많은 것들을 체크할 수 있다. ximage_lbfgs.o:ximage_lbfgs.f90 timpose.o $(cmpl) $(Ocft) ximage_lbfgs.f90 ximage_lbfgs.o를 만들기 위해서는 timpose.o가 먼저 필요하게 된다는 것입니다. 또는 timpose.f90이 수정되었 다면 다시 컴파일되어야 합니다. 이러한 변화에 발 맞추어서 ximage_lbfgs.o도 반드시 다시 컴파일되어야 한 다는 의미입니다. 여기에서, 재미있는 것을 얻을 수 있습니다. 프로그램의 버전이 올라갈 때마다 (또는 프로 그램 개발 중에), 변화가 상당히 국소적으로 일어남을 알 수 있습니다. 즉, 프로그램의 전체적인 변화가 분명 히 있음에도 불구하고 구체적인 소스코드 상에서의 변화는 상당히 국소적인 변화로 종결될 수 있음을 보여 주고 있습니다. 이것이 바로 대상중심 프로그램의 장점 중의 하나입니다. 즉, 전체 설계도를 보지 않아도 된 다는 것입니다. 프로그램이 크게 바뀔 때마다 모듈의 형태로 추가될 수 있다면 이 또한 보기 좋을 것입니다. 한 가지 재미있는 것은 아래와 같다. 여러 모듈 속의 변수들을 이제는 너무 쉽게 접근하고 바꿀 수 있어지니 까 '모듈 사용의 도'를 넘어 오히려 전체 프로그램을 못 읽어 버리게 만들 수 있다는 사실이다. 찬찬히 생각하 면서 모듈의 기능과 전체 프로그램의 흐름을 다시 생각해 볼 필요가 있다. 이러한 경우 새로운 모듈들의 재정 립이 필요하다. IMPLICIT NONE에 관해서: 어떠한 것도 묵시적이지 않다고 선언하는 것이다. 즉, 모든 변수들이 구체적인 형을 선언해 주어야 한다는 것이다. 이렇게 선언한 이상. REAL*8, INTEGER, CHARACTER*3, LOGICAL, COMPLEX*16 과 같이 모든 변수에 그 해당 변수형을 선언해야 한다. 이러한 것은 프로그램 작성시 발생한 TYPO들을 잡아주는 기능을 한다. 조금 귀찮지만 안전한 코드 생성에 일조한다. 더욱이, 모듈 개념을 사용할 때는 꼭 있어야 한다는 생각이 든다. 왜냐하면, USE를 사용하는 경우, 사용하는 모듈 속에 변수들이 존재하 는지를 컴파일하는 단계에서 강제로 확인 시켜준다. 이전에 많이 사용했던 IMPLICIT REAL*8 (A-H,O-Z) 즉, 시작하는 변수들이 A-H 또는 O-Z로 시작하면 실수 형임을 선언하는 어찌 보면 너무나도 광범위한 선언은 프로그램 작성과 관리차원에서 그렇게 좋은 선택은 아니다. 복잡한 프로그램을 만들다 보면, 항상 사용하는 변수들이 선언되게 하는 implicit none의 힘을 느낀 다. 뻔한 내용의 실수들을 줄일 수 있도록 도와주는 힘. 뻔한 실수들을 잡아내기가 더욱 어려울 때가 종종있 다. 전혀 예측하지 않은 형태의 착오이기 때문이다. 변수 선언은 INTEGER, DIMENSION(:), ALLOCATABLE, SAVE, PRIVATE :: ivari_name 처럼 사용될 수 있다. 너무 길어 (5 of 114) 오전 10:21:34

6 지는 면이 있지만, 변수하나의 특성을 표현하기 위해선 어쩔 수 없다. 실제로는 이와같이 하지 않는 경우가 많다. INTEGER, ALLOCATABLE :: ivari_name(:) 정도로 해두고 private와 save는 보통, 상당히 앞 쪽에서 미 리 따로 선언하는 경우가 많다. 또한 재 선언을 통해서 public 또는 private를 확인해 줄 수 있다. implicit none private save integer... real*8... logical... public ::...(부)함수이름, 변수이름, (가변)배열처럼 되는 경우가 많이 있다. USE를 사용할 때에도 USE xmodulex, ONLY: ivari,jvari,isubp,ifunc 처럼 ONLY를 사용해서 구체적으로 함수와 변수들을 명시해주 는 것이 보기에 좋고 프로그래밍에서 실수할 확률을 줄일 수 있도록 도와 준다. 프로그램 시작과 함께 변수값 이 할당될 수도 있다. 이 경우는 integer :: iiijjj=10, kkklll=100 logical :: on=.true., off=.false. real*8 :: qqq=1.10d0, ttt=2.20d0 real*8 :: aaa(-1:1)=(/ -1.d0, 0.0d0, 1.d0 /) 통상의 구조는 어떠한가? (제법 많을 일들을 수행하는 경우) = (특정 데이터의 크기를 모듈 외부와의 통신을 통하여 정하고, 필요하면 디스크로부터 정보를 읽어 들이고, 또는 생성하고, 연관된 계산들을 수행하고, 주 요 데이터들을 모듈 외부로 보내고/받고, 특정 배열을 소멸시키고, 주요 정보를 디스크에 저장하고, 프린트하 고,...) module common_data implicit none private save integer iii1,jjj1 real*, allocatable :: aa1(:,:),bb1(:) logical first_call 첫 번째로 모듈이 사용되는 경우가 중요한 경우가 있다. 이를 위해서 전역변수로 선언해 둘 수 있다. 초기화를.true.로 해두고, 나중에.false.로 바꿀 수 있다. character*2, allocatable :: ccc1(:)... public :: aa1,bb1,iii1 end module common_data module common_data2 implicit none private save... public ::... end module common_data2 program abcd implicit none (6 of 114) 오전 10:21:34

7 contains subroutine aaa USE common_data implicit none... end subroutine aaa subroutine bbb USE common_data2, ONLY : w,y,z implicit none... end subroutine bbb subroutine ccc USE common_data implicit none... end subroutine ccc subroutine ddd USE common_data2, ONLY : x,y,z implicit none... end subroutine ddd subroutine eee implicit none... end subroutine eee function fff implicit none... end function fff end program abcd implicit none 뒤에는 각종 변수, (가변) 배열들이 선언될 수 있다. 상당히 짧은 프로그램의 경우 위와 같은 형 식으로 해결될 수 있다. 일반으로, 사용되는 모듈들의 숫자가 많아 질수록 대부분의 경우 여러 개의 파일들 을 사용한다. 보통 한 파일 속에 모듈 하나를 넣어 두는 경우가 많다. 위의 경우에는 모듈들은 데이터 저장 창 고로 사용된 예이다. 하지만, 모듈은 특정 데이터들을 보관할 뿐만 아니라 관련된 데이터 처리를 할 수 있도 록 설계되었다. 즉, 모듈= {데이터들}+{데이터 처리 부함수들}이기 때문이다. 모듈이 변수들 없이 함수들만 포함(contain)할 수도 있다. 반대로 변수들만 가질 수도 있다. program ssss USE xxx, ONLY : aa1,bb1,pbc USE yyy, ONLY : qq,pp,ssq USE www implicit none (7 of 114) 오전 10:21:34

8 REAL*8 ss REAL*4 t_s INTEGER itemp,itemq,irate CHARACTER*8 fnnd ; CHARACTER*10 fnnt real*8, allocatable :: vector(:) integer ierr CALL DATE_AND_TIME(date=fnnd,time=fnnt) write(6,'(1x,a10,2x,a8,2x,a10)') 'date,time ', fnnd,fnnt CALL SYSTEM_CLOCK(itemp,irate) write(6,*) epsilon(ss),' epsilon for real*8' write(6,*) epsilon(t_s),' epsilon for real*4' E-016 epsilon for real* E-07 epsilon for real*4 write(6,*) floor(4.8d0) write(6,*) floor(-4.8d0) 4-5 write(6,*) huge(ss) write(6,*) huge(t_s) E E+38 write(6,*) minexponent(ss) write(6,*) minexponent(t_s) write(6,*) nint( ) write(6,*) nint( ) -7-6 가장 큰 수를 찾고자 할 때 max0(j,k,l,...) amax0(a,b,c,d...) ifix(aa)버림을 실행하여 정수로 반환 min0(i,j,k,...) amin0(a,b,c,...) sign(x,y) 는 abs(x)*(y 값의 부호)를 반환한다. mod() MOD (3.0, 2.0) (8 of 114) 오전 10:21:34

9 MOD (8, 5) 3. MOD (-8, 5) -3. MOD (8, -5) 3. MOD (-8, -5) -3. int() 정수로 변환 INT (-3.7) -3. nint() Nearest integer. NINT (2.789) --> 3. NINT (2.123) --> 2. aint() (n/2)*2 n=11 8 n=11 10 n=-9-8 포트란90에서의 수학 DOT_PRODUCT(vector1,vector2) MATMUL(matrix1,matrix2) MAXVAL() MINVAL() SUM() ==.eq. /=.ne. 비교: C 언어, 파이썬에서는 =을 사용한다. >=.ge. <=.le. 위의 경우처럼 '='이 관련된 경우에서 =이 항상 뒤쪽에 위치한다. >.gt. <.lt. 논리연산에는.not.,.and.,.or.가 있다..not. ( (la == lb).and. (.not.(lc /= ld) ) ) 처럼 논리 연산에 직접 사용될 수 있다. if(.not. allocated(vector))then allocate(vector(100000), stat=ierr) if(ierr /= 0) then write(6,*) 'problem' stop endif endif allocate(vector(-89:100)) 처럼 정의될 수 있다. -89, -88, -87,...98, 99, 100. 즉, allocate(vector(100)) 는 크기가 100인 vector가 생성됨을 의미한다. 물론, 1부터 시작하여 100까지 간다는 말이다. allocate(vector(1:100))과 동 (9 of 114) 오전 10:21:34

10 등한 표현이다. 디폴트가 1부터 시작한다는 이야기입니다.... if(allocated(vector)) deallocate(vector) CALL SYSTEM_CLOCK(itemq) write(6,'(2e15.4,2x,a9)') float(itemq-itemp)/float(irate)/60.,float(itemq-itemp)/float(irate)/3600.,' min or h' end program ssss module xxx implicit none private save 모듈 내에서 공통으로 사용되고 저장될 변수들...global object 선언 (전역 객체 선언) 이곳에, 모듈 내에서 정의되는 함수 이름(real*8, integer, logical 등을 선언할 필요는 없다. 모듈 내에서 각자 정의될 때만 필요하다. logocal first_call... public ::...변수, 가변 배열, 부함수 (subroutine, function) 이름들... contains subroutine xxx_initialize(...) implicit none first_call=.true....지역 변수들 local variables end subroutine xxx_finalize subroutine xxx_finalize(...) implicit none...지역 변수들 end subroutine xxx_finalize subroutine xxx_do_something(...) USE xxx, ONLY : xxx1,yyy1,sub_name,fn_name implicit none...지역 변수들 주의할 점: 지역변수들을 선언할 때, 전역변수와 이름이 같은 것을 사용해 버리면 어떻게 될까? 답은 간단하 다. 이름이 같은 것조차 의미없는, 진정한 의미로, 한번 지역 변수는 그 지역에서 지역변수이다. 하지만, 이렇 게 복잡하게 프로그래밍할 필요가 없다. 혼돈을 피하기 위해서, 지역변수들은 될 수 있는 한 전역변수와 같 은 이름을 사용하지 않도록 하는 것이 바람직하다. 모듈 속의 서브루틴에서 이 모듈 속에 속하지 않는 외부의 서브루틴, 함수들도 당연히 사용할 수 있다. 독립 적으로 선언된 모듈 속의 내용(변수들, 부함수들)이라면 당연히 USE를 사용해야 한다. 예전에 만들어진 포트 란77 프로그램과 함께 연동될 경우, 특히, common문이 있는 경우에도 common문을 그냥 사용해도 된다. 왜냐 하면 포트란90은 포트란77을 다 이해하기 때문이다. end subroutine xxx_do_something subroutine xxx_print USE www USE zzz, ONLY : sfa 모듈 속의 부함수들은 인수를 가지지 않을 수도 있다. 또 다른 모듈의 변수,함수를 사용할 수 있다. 또 다른 모듈의 변수,함수를 사용할 수 있다. (10 of 114) 오전 10:21:34

11 implicit none... end subroutine xxx_print function x_test() implicit none real*8 x_test 모듈 속의 함수 선언 :이곳에서 함수 선언하는 것만으로 충분하다. 다만, 전역 변수 선언 공간에 서, 이 친구가 private, public이라는 것만 밝혀두면 된다. end function x_test end module xxx 마찬가지로 module yyy도 고유의 모듈 형식을 가질 수 있다. module www, zzz등도 있어야 컴파일이 된다. 모듈 내부에서만 저장되거나 생성되는 변수, 내부에서만 이용되는 부함수 모듈 외부에서 (USE 명령을 통해서 사용) 참조가능한 변수, 부함수 등에 유의하여 작성한다. MODULE x_mo_du_le IMPLICIT NONE PRIVATE 디폴트를 어떻게 잡느냐의 문제; 기본이 private이고 public에 해당하는 것들은 표식적으로 명시해야 함을 나타내고 있다. SAVE INTEGER :: ia,ib,ic,id,ie,if... REAL*8, ALLOCATABLE :: arr(:,:), brr(:,:), crr(:,:)...(특정한 순간에 크기가 할당될 배열들이다.)(계산/저장 후) (특정 순간에 할당된 정보가 없어질 수 있다.) PUBLIC :: abc_initial, abc_final,variable_1,array_1,ia,ib,...필요한 것만 직접 적어둔다. 이 영역에서는 공통으로 사용되는 변수,배열,가변배열들을 선언해준다. 아래에 선언될 함수, 서브루틴들에 의해서 참 조될 뿐만 아니라 변화될 수 있다. 물론, 외부에서 참조되고 변화될 수 있다. 외부에서 사용하는 경우 USE라는 명령을 사용하여 그곳에서 특정 모듈 내용을 사용하게 된다. 여기서 주목할 것은 private 과 public이다. 처음에 private으로 시 작했기 때문에 이 영역에서 선언된 것들은 이 모듈이 사적으로 사용한다는 뜻이다. 즉, 이 모듈 밖에서는 이 변수가 보 이질 않는다. 따라서 예기치 않은 상황의 연출로부터 바뀔 염려가 없다. 완전히 보전된다. 다만, 뒤에 마음이 바뀌어 서 구체적으로 public (즉, 밖에서 사용할 수 있는 subroutine 이름이나 배열, 정수, 실수 등을 )을 이용하여 선언해준다. 즉, 기본 사적사용에 공적사용 명시인 것이다. 위에서 선언된 save는 기본적으로 모든 변수들이 이 모듈이 사용되는 동안 그 값이 계속 저장된다는 것이다. 연이어서 이 모듈을 사용할 때 이전 값들이 그대로 유지되어야 한다고 선언하 는 것이다. private, public을 사용하지 않으면 public으로 간주된다. 따라서 첫 줄에 private를 사용하지 않았다면,마지막 에 private를 이용하여 현 모듈에 있어서 사적인 것들의 리스트를 만들어 두면 좋다. 예를 들면, 아래와 같다. public :: aa, bb, cc public :: aa_subroutine, bb_subroutine, cc_function 이렇게 public으로 선언된 것들은 원하는 모듈 또는 부 함수들에서 USE, ONLY들을 이용하여 선택적으로 사용할 수 있다. CONTAINS 뒤에 따르는 함수들을 나열하기 전에...;물론, 서브루틴, 함수들은 위에서 선언된 변수들을 중점적으 로 관리할 서브프로그램들일 것이다. 즉, 특정 대상 중심의 서비스. 국소 변수들로 argument형식으로 외부와 통신할 수 있다. 하지만 주요 정보들은 이곳에서 계산되고, 이곳에 저장되면 될 것이다. 일반으로, 본 모듈 내에서만 정의되고, 사용되는 변수들과 부함수들도 있다. (이들을 private이라 고... ) 모듈 외부에서 접근/사용 가능한 변수들과 부함수들===>public이라고... 외부와 통신하는 부함수들은 자체적으로 인수들을 가질 수 있다. 국소적인 인수들이다. 위에서 선언된 모듈 내의 변수들과는 명목적으로 상관이 없는 변수들을 인수로 가지고 있다. 인수들을 가지지 않는 부함수들도 있다. 이 경우 이미 주요한 변수들이 위에서와 같이 정의된 경우, 원하는계산만 수행할 경우 특별 (11 of 114) 오전 10:21:34

12 히 인수가 없을 수도 있다. 통상 외부 함수 (주 프로그램,서브루틴, 함수에서 접근 가능한 subroutine들이 예를 들어 3개 정도는 있을 수 있다. 물 론, 충분히 많은 수 일 수 있다.; 모듈 내부에서 만 사용되는 변수/부함수들도 있을 수 있다. 필요한 일들을 충분히 분업 화 시킬 수 있다는 의미이다. 위에 선언된 모듈 내에서만 정의되어 보호되는 변수들을 전문으로 취급하는 부 함수들 을 많이 만들 수 있다. 특정한 데이터들만을 전문으로 다루기 때문에 이러한 부 프로그램들은 통상 private형식인 것 이 자연스러운 일일 것이다. 이와는 반대로, 외부와 직접 정보를 교환하는 부프로그램은 당연히 public 형식일 것이고 (다시 말해서, 외부에서 접근가능한 부 프로그램임을 의미한다.) 필요에 따라서, 용도에 따라서, 인수들 (변수, 배열을 가질 수 있다. SUBROUTINE abc_initial(n_problem_size,...) IMPLICIT NONE 통상 main함수로부터 불려지면서 (따라서 public이어야 한다.) 현 모듈과 관련하여 필요한 메모 리를 확보한다. 또는 변수 초기화를 시도한다. INTEGER n_problem_size 일반적으로 주 (main) 프로그램에서 인수들을 통해서 문제의 크기를 정해줄 수 있 다. 직접 필요한 정보를 읽어들여도 됨. INTEGER i1_local,i2_local REAL(8) r1_local,r2_local... 일반적으로 총체적인 입력들을 정리해 줄 수 있다. 실질적으로 풀려는 문제가 구체화되고, 정의된다. 사실, 관련 모듈 의 특성화가 시작된다. ALLOCATE(aa(n_problem_size)) ALLOCATE(bb(n_problem_size)) 모듈 내에서 언제든지 접근가능하다. 물론, public이면 외부에서도 접근가능 함.... END SUBROUTINE abc_initial SUBROUTINE abc_final() IMPLICIT NONE INTEGER i1_local... main함수로부터 불려지면서 (따라서 public이어야 한다.) 불 필요한 메모리를 지운다. REAL*8 r1_local,r2_local... 일반적으로 이 모듈과 관련된 일들을 종료하기 전 총체적인 출력 정리를 해 줄 수 있다. 필요한 정보들 출력, 예를 들 면. DEALLOCATE(aaa) DEALLOCATE(bbb)... END SUBROUTINE abc_final SUBROUTINE aaaaa1() 서브루틴, 함수의 인수가 하나도 없을 경우도 있을 수 있다. IMPLICIT NONE INTEGER i1_local,i2_local REAL*8 r1_local,r2_local...구체적으로 이 모듈 내에서 할 일 들과 연관이 있는 일을 최대한 분업화할 수 있도록 한다. 여러 개의 서브루틴과 함수를 정의한다. 모듈 내에서는 많은 인수들이 필요하지 않다. 왜냐하면 모듈 내에서 대부분 의 변수들이 접근가능하기때문이다. 인수가 하나도 없는 함수도 물론 가능하다. 즉, subroutine unaaa, 부를 때 call unaaa 형식이 가능하다. argument를 이용 하여 외부와 최소한의 통신은 가능하도록 하는 것이 일반적이다. 모듈 내에서 사용되는 다양한 변수들은 외부 모듈 들 속의 변수들 또는 다른 모듈들과의 접촉을 차단한 채 그들만을 위하여 특별히 만들어진 부 프로그램들에 의해서 (12 of 114) 오전 10:21:34

13 다루어진다. 모듈의 개념을 도입하면 앞서 이야기한 특정 대상을 중심으로 프로그램들이 쉽게 만들어진다. 상황이 이 러하니 인수들을 많이 가지지 않는 것이 당연하다. END SUBROUTINE aaaaa1 SUBROUTINE internal_use1 IMPLICIT NONE INTEGER i1_local,i2_local REAL*8 r1_local,r2_local 현 모듈 내에서만 사용되는 부함수가 일반으로 많이 있을 수 있다. 인수가 있을 수도 없을 수도 있다. 물론, 목 적과 설계에 따라서 달라질 수 있다. 모듈 속에서 선언된 전역변수들만 집중적으로 취급하는 경우는 모듈 속 서브루틴, 함수들은 인수의 숫자가 적어질 수밖에 없다. 물론, 모듈 외부와 통하는 경우에는 일반으로 인수들 의 숫자가 증가할 가능성이 많을 것이다. END SUBROUTINE internal_use1 FUNCTION qqq() 서브루틴, 함수의 인수가 하나도 없을 경우도 있을 수 있다. IMPLICIT NONE REAL*8 qqq INTEGER iq_local, REAL(8) qx_local... RETURN 없어도 된다. subroutine의 경우에도 return없어도 된다. END FUNCTION qqq 이 영역은 필요한 함수나 서브루틴들을 만들어서 두는 곳이다. 통상의 인수 (arguments)들은 최소화되어 몇 개밖에 되 질 않을 것이다. 왜냐하면 이 모듈 내에서 공통으로 사용되면서 계속 저장되는 (save) 변수들은 이미 모듈 전역 변수들 로 정의되어 있기 때문이다. 또한 필요에 따라서 데이터들을 저장하는 배열들을 생성시킬 수 있다. 이 경우에는, 모듈 내의 함수들은 외부에서 집중적으로 불려지는 것(주요 함수)과 한 번 정도 불려지는 것들 (_initial, _final) 등으로 분류될 수 있다. 또한 현 모듈 내에서만 불려지면서 사용되는 서브 프로그램들이 있 다. 모듈 속에 선언된 전역 변수들도 마찬가지이다. 모듈 내에서만 사용되는 변수들도 일반적으로 있을 수 있다. END MODULE x_mo_du_le 이 모듈을 사용하려면 사용하고자 하는 함수, 주 프로그램, 또는 서브루틴 내에 USE x_mo_du_le (물론, USE x_mo_du_le, ONLY : xxx 처럼 선언하는 것이 더 강력한 프로그래밍 습관이다.)이라고 선언해 준다. 하지만, 그 선언되는 위치는 주의해야 한다. 보통 implicit none 위에서 선언한다. 엄청나게 앞에서 선언되는 것이다. 현 함수에서 선언되는 변수들 보다 항상 앞선다는 것이다. 여러 개의 모듈들은 동시에 사용할 경우 연속해서 한 줄 단위로 추가하여 적어두면 될 것이다. 이렇게 하면 앞서 이야기한 public이 허용하는 부함수들 또는 (가변) 변수들이 현 함수 내에서 접근/조작이 가능해지는 것이다. 물론, 각 모듈들 내에서 private변수들은 현 함수에서 볼 수가 없다. 만약, 현 함수에서 implicit none을 사용했다면, 프로그램들을 컴파일할 때 즉시, 잘못 사용되고 있는 변수들을 찾아낼 확률이 상당히 높아지는 것이다. (물론, 이것이 프로그래밍에서 함정이 될 수 도 있다. 우연히 지역변수를 선언하고 사용하고 있는데, 이것이 다른 모듈 속에서 온 것으로 착각하고 프 로그램을 하는 경우이다.) 이러한 자체 체크 기능은 상당히 중요한 의미를 부여한다. 여러 명이 동시에 일을 하다보면 뜻하지 않은 일들이 많이 발생하게 된다. 혼자 프로그램을 만들어도 마찬가지이다. 항상 그 내용과 이름을 다 기억하고 있을 수는 없을 것이다. 시간이 지나고 다음에 다시 조금 수정할 때에도 마찬가지이다. 안전사고의 미연 방지는 아무리 강조해도 지나치지 않다. 포트란 프로그램에서 `시스템 콜하기`; 파이썬과의 결합?포트란 프로그램에서 CALL SYSTEM('application. (13 of 114) 오전 10:21:34

14 py')처럼 파이썬 스크립트를 직접 불러서 사용할 수도 있다. 이러한 경우 call system이전에 필요한 자료들을 파일로 정리해둔 다음 시스템으로 잠시 나간 다음 시스템에서 준비된 파이썬 스크립트가 자료를 사용할 수 있게 설계할 수 있다. 물론 외부의 프로그램들 (스크립트에서 준비) 등으로 자료처리 또는 계산을 수행한 다 음 다시 포트란 프로그램으로 되돌아 오는 방법이 있다. 물론 돌아오기 전에 포트란 프로그램에서 자료를 읽 어들일 수 있도록 파일들을 미리 만들어주면 좋겠다, 물론 파이썬 스크립트가 관리를 해야 할 것이다. 마지막 으로 포트란 프로그램에서 준비된 자료들을 읽어들여서 계속하여 포트란 프로그램으로 계산을 수행할 수 있 다. 자료 적어두기: 파일열기/적기 ==포트란 명령으로 수행 call system('python_driver.py') ==포트란에서 빠져 나와서 파이썬 스크립트 실행하기. 자료 읽어들이기: 파일 열기/읽기 ==포트란 명령으로 수행 유닉스 프롬프트 상에서 아래와 같이 abc.exe라는 실행화일을 실행시키면 nohup abc.exe >output_file_name & 사용자 계정에서 로그아웃해도 프로그램 실행이 죽지 않는다. derived-type statement type sphere real*8 :: cx,cy,cz,radius end type sphere type (sphere) :: a_ball, b_ball 인터넷에 공개된 포트란90 관련 자료들 lapack, blas, atlas (14 of 114) 오전 10:21:34

15 간단한 모형들 프로그램 분석을 위해서 프로그램을 프린터할 때 a2ps -o output.ps <input.f90 처럼 a2ps프로그램을 이용하면 보기 좋은 PS파일이 생긴다. 모형 (1) 하나의 간단한 변수 저장 창고, globy라고 이름 붙여진 정보 저장창고이다. 접근 가능한 모든 부 프 로그램들에 의해서 참조 및 변수값의 변화가 가능하다. save로 저장되기 때문에 계속적으로 변화가 업데이 트 된다. 이러한 데이터 창고가 여러 개 있을 수도 있다. 이 정도는 FORTRAN77의 named common 정도로 해 석될 수 있다. {common /aaa/ ig1,ig2,ag1,ag2,ag3} 물론, 크기가 정해진 배열이 포함될 수도 있다. 배열의 크기 가 미리 정해지지 않은 경우는 아래에 있는 예제들에서 다룬다. 실행할 때 그 크기를 정하고 계속사용한다. MODULE globy implicit none real*8, save :: global_variable1,global_variable2 이러한 경우 default를 따라서 모두 다 public하다. integer, save :: iglobal1,iglobal2 END MODULE globy 하나의 데이터 저장 창고이다. 가장 간단한 예이다. 필요에 따라서 아래의 예에서와 같이 각종 형태의 가변 배열도 선언할 수 있다. integer, allocatable :: ijk(:) MODULE globy implicit none PUBLIC SAVE real*8 global_variable1,global_variable2 integer iglobal1,iglobal2 END MODULE globy PROGRAM main1 USE globy implicit none보다 위에 위치한다. 컴파일할 때도 마찬가지이다. "globy"를 가지고 있는 상황에서 컴파일 이 진행될 수 있다. USE globy, ONLY : iglobal1 IMPLICIT NONE REAL*8 a1,a2,a3 INTEGER ia1,ia2,ia3 READ(5,*) igloval1, iglobal2 READ(5,*) global_variable1,global_variable2...중략... END PROGRAM main1 위와 같은 형식의 여러 가지 모듈들을 사용하여 필요한 프로그램 단위들에 걸쳐서 사용할 수 있다. 즉, 필요 한 프로그램 단위에 USE를 사용하여 데이터를 공유하고 변화시킨다. (15 of 114) 오전 10:21:34

16 MODULE v_globy arrays to be shared... IMPLICIT NONE PUBLIC SAVE REAL(8), ALLOCATABLE :: v_global(:,:) INTEGER, ALLOCATABLE :: i_global(:) INTEGER ii,jj,kk REAL(8) aa,bb,cc, ddd(100,10) END MODULE v_globy 가변 배열을 포함하는 형식이다. 포트란77과 비교하면 상당한 수준의 발전을 느낀 다. 프로그래머가 필요한 메모리 할당량을 구체적 프로그램 실행 시에 조절할 수 있게 해 준다. 더 이상 필요 하지 않을 경우 소멸시킬 수 있다. 사용하는 입장에서는 가변 배열의 길이가 할당되었는지를 확인하고 사용 하면 된다. 예를 들어, 아래와 같이. if(.not. allocated(v_global)) then allocate(v_global(ii,jj)) endif 마찬가지로 필요가 없을 때 해당 모듈을 사용하는 적절한 프로그램 단위에서 deallocate(v_global)할 수 있다. 물론, 할당된 배열인지 구분을 위해서 allocated()함수를 사용할 수 있다. 여기에 제시된 모듈의 형식은 아주 단편적인 것이다. 보다 전문적으로 사용되는 모듈의 형태를 알아 둘 필요가 있다. 할당되었는지를 확인하고 또한 필요없을 경우 할당량을 없애 버리는 작업: if(allocated(v_globy)) deallocate(v_globy) 모형 (2) MODULE xmodule IMPLICIT NONE PUBLIC 이번엔 기본 public입니다. SAVE REAL*8, allocatable :: ak(:,:,:),qq(:,:,:),gradix(:,:,:) "::" 가 반드시 있어야 한다. REAL*8, allocatable :: force(:,:,:),vofqj(:),xmass1(:) REAL*8, allocatable :: eksq(:,:,:),fk(:,:,:) REAL*8 :: tau,xmu,xmass,onma,object,etarget,detar,amp,fdeltak,xnu,ttarget REAL*8 penalty0,penalty1,penalty2,penalty3 CHARACTER*2, allocatable :: isymbol(:) INTEGER np,nk,natom,itang,iprint,itmax_relax INTEGER :: iatom_1,iatom_2,iatom_3 INTEGER iatom_1 처럼 선언해도 된다. 즉, :: 없이 선언할 수 있다. LOGICAL lfirst_xim,llast_xim LOGICAL :: lneb,lcho,lpar,lclimb,lkntar LOGICAL :: lperpath REAL*8 :: aperc,pamp PRIVATE :: qq,gradix,force,eksq,fk 기본 public인데, 이들은 private임을 명시하고 있다. 이들은 현재의 모듈 내에서만 사용될 수 있다는 뜻이다. 다른 곳에서는 보이지 않는다. 부 프로그램 이름, 변수 이름등이 private라는 특성 을 가질 수 있다. CONTAINS SUBROUTINE alphaomega(nvar,xarray) USE impose_mod 또 다른 모듈을 사용할 때, 제일 앞 쪽에서 선언한다. implicit none보다 더 앞서서... IMPLICIT NONE INTEGER nvar REAL*8 xarray(nvar) REAL*8 pi,tmv(3),uu,us,uf,du,del REAL*8 tmp,ej,vofqj_max,vofqj_min REAL*8, ALLOCATABLE :: work1(:),work2(:) INTEGER, ALLOCATABLE :: ind1(:),ind2(:) REAL*8 rmsvalue REAL*8, ALLOCATABLE :: xtmp(:),ytmp(:),ztmp(:) REAL*8, ALLOCATABLE :: xtmq(:),ytmq(:),ztmq(:) (16 of 114) 오전 10:21:34

17 INTEGER k,j,kp,iatom,natomold,nkold LOGICAL lexist CHARACTER*2 ichar,jchar LOGICAL latomic_a,latomic_b REAL ranmar...중략... 첫 번째 call인지 마지막 call인지를 판단하여 각각 allocate,deallocation을 실행하고 입력/출력 정리를 행한다. RETURN END SUBROUTINE alphaomega SUBROUTINE genthetax(theta,posi,grad,nn) IMPLICIT NONE INTEGER nn REAL*8 theta REAL*8 posi(nn),grad(nn) 외부와 직접 통신하는 변수들입니다. 명목적으로는, 형식상 모듈 내의 변수들과 전혀 상관 없습니다. REAL*8, allocatable :: works(:) REAL*8, allocatable :: workt(:) REAL*8 tmv(3),tmw(3),tang(3),tnu(3),tangp(3),tangm(3) 3-d problem REAL*8 cc,ac,uu,us,uf,du REAL*8 tmr,tmq,ej,pi,del,tmp,fphi,e_ref,e_max INTEGER iatom,j,k,kp REAL*8 tjaver,zz 이들은 국소 변수들 (local variables)이다. 모듈 내의 다른 함수들/전역변수들과는 전혀 상관없다. 설 사 이름이 같다고 하더라도 전혀 상관없다. decoding ---[ DO iatom=1,natom kp=3*nk*(iatom-1) DO k=1,nk ak(iatom,1,k)=posi(kp+3*(k-1)+1) ak(iatom,2,k)=posi(kp+3*(k-1)+2) ak(iatom,3,k)=posi(kp+3*(k-1)+3) ENDDO ENDDO 중략 posi,grad는 주 함수에서만 정의된 일차원 변수이다. 사실 상 이 변수들이 현 모듈에서는 2차원 배열등이 될 수도 있 다. 예를 들면, 동적인 변수 최적화의 경우 lbfgs같은 루틴들은 1차원 배열을 입출력으로 요구하는 반면, 현 모듈 내에 서는 이것이 원자의 위치를 나타내는 변수로서 2차원 배열일 수 있다. 이럴 경우 적절한 이들의 관계를 잘 정의하여 사용할 수 있다. decoding과 encoding을 본 루틴의 맨 앞 부분과 맨 뒤 부분에 정의해 두면 아주 편리하다. END SUBROUTINE genthetax 중략 (17 of 114) 오전 10:21:34

18 END MODULE xmodule 하나의 전형적인 모듈 모형을 아래에 표시했다. module xxx implicit none private save... public :: yy,zz,...얼로케이터블 배열명, 변수명, 함수명, 부함수명,...: 이 모듈 밖에서 USE를 사용해서 사 용할 수 있고, 변화시킬 수 있는 것들. contains subroutine xxx_initial(...) implicit none 지역 변수 선언 (가변배열, 배열, 실수, 정수,...) 이곳에서 독립적으로 선언문을 만들 수 있다. 선언문 형식으로 배열, 가변 배열, 실수, 정수등을 선언할 수 있 다. 따라서, 이 들은 이 부 프로그램에서만 사용되는 지역변수들이다. 예를 들면, 특정한 인수들을 이용하여 현 모듈에서 필요로 하는 객체들을 준비한다. allocate(...,...)... end subroutine xxx_initial subroutine xxx_final implicit none 지역 변수 선언 (가변배열, 배열, 실수, 정수,...) 모듈을 종료하는 프로그램 deallocate(..,..)... end subroutine xxx_final subroutine xxx_maijor implicit none 지역 변수 선언 (가변배열, 배열, 실수, 정수,...) 현 모듈이 추구하는 주된 계산들을 행한다.... end subroutine xxx_major subroutine xxx_subs1 현 모듈 밖에서 결코 사용되질 않는다면, 모듈에서, 전역변수 선언하는 곳에서, private라고 명시해줄 필요가 있다. implicit none 지역 변수 선언 (가변배열, 배열, 실수, 정수,...) 몇 가지의 보조 프로그램들이 필요할 수 있다.... end subroutine xxx_sub1 (18 of 114) 오전 10:21:34

19 subroutine xxx_subs2 implicit none 지역 변수 선언 (가변배열, 배열, 실수, 정수,...) 현 모듈 내에서만 사용될 수도 있다. 모듈 속에서 공통으로 사용가능한 정보들은 그냥 접근할 수 있다. 변화 시킬 수도 있다. 전역변수이기 때문에.... end subroutine xxx_sub2 function ss 현 모듈 밖에서 결코 사용되질 않는다면, 모듈에서, 전역변수 선언하는 곳에서, private라고 명 시해줄 필요가 있다. implicit none 지역 변수 선언 (가변배열, 배열, 실수, 정수,...)... end function ss end module xxx 모형 (3) MODULE example1 IMPLICIT NONE PRIVATE 기본적으로 사적인 변수들만 있다. 예외는 아래에서와 같이 public으로 확실히 명시한다. SAVE 아래에 표시될 모든 변수들은 본 모듈 속에 안전하게 저장될 것이다. PUBLIC :: ex_action, ex_initial,ex_final INTEGER :: ia,ib,ic,ions 중략 FUNCTION dist(v1,v2) REAL*8 :: V1(3,ions),V2(3,ions),V12(3,ions),dist V12=V2-V1 ; dist=sqrt(sum(v12**2)) END FUNCTION dist FUNCTION unit(v1) REAL*8 :: V1(3,ions) REAL*8, dimension(3,ions) :: unit unit=v1*(1.0d0/sqrt(sum(v1*v1))) END FUNCTION unit 중략 모형 (4) lbfgs 알고리즘을 불러서 사용자 정의 함수 최소화하기: 사용자가 정의한 모듈 속 에서 사용자 정의 서브루틴 (genthetax)과 같이 아래에 있는 서브루틴이 나란히 정의되면 된다. (함수를 최소 화할 때 lbfgs는 상당히 유용한 루틴이다. 현실적으로 가장 빨리 함수를 최적화하는 루틴 중 하나로 보인다. SUBROUTINE lbfgs_header(theta,posi,grad,nn) written by In-Ho Lee, KRISS, Jan (19 of 114) 오전 10:21:34

20 IMPLICIT NONE INTEGER nn REAL*8 theta,posi(nn),grad(nn) integer n_local,m,msave real*8, allocatable :: diag(:),w(:) real*8 eps,xtol,gtol,t1,t2,stpmin,stpmax integer iprint(2),iflag,icall,mp,lp,j,nwork logical diagco The driver for LBFGS must always declare LB2 as EXTERNAL external lb2 common /lb3/mp,lp,gtol,stpmin,stpmax LBFGS는 포트란77로 만들어진 프로그램입니다. 바꾸고 싶은 생각이 전 혀 없습니다. 그냥 사용하는 것이 목표. 함수값과 gradient를 이용해서 함수 최소화시키는 루틴입니다. quasi Newton 방법으로 알려져 있습니다. n_local=nn m=5 ; msave=7 ; nwork=n_local*(2*msave+1)+2*msave allocate(diag(n_local),w(nwork)) iprint(1)= 1 ; iprint(2)= 0 We do not wish to provide the diagonal matrices Hk0, an therefore set DIAGCO to FALSE. diagco=.false. ; eps= 1.0d-3 ; xtol= 1.0d-16 ; icall=0 ; iflag=0 eps 가 작을수록 더욱더 정확한 함수 최소화를 실행한다. 20 continue call genthetax(theta,posi,grad,nn) 함수값과 gradient를 계속 계산하도록 되어있습니다. 함수 최소화와 관련된 작 업은 20 번 루프와 lbfgs루틴이 알아서 수행함. write(6,'(4e14.7,1x,e10.4,2x,e11.5)') theta,penalty0,penalty1,penalty2,onma,etarget grad의 부호가 -grad로 되면 계산이 안된다. 즉, lbfgs루틴은 정확히 gradient를 요구한다. -gradient, 즉, force를 사용하 면 안된다. call lbfgs(n_local,m,posi,theta,grad,diagco,diag,iprint,eps,xtol,w,iflag) if(iflag.le.0) go to 50 icall=icall + 1 We allow at most 2000 evaluations of Function and Gradient if(icall > 2000) go to 50 go to continue deallocate(diag,w) END SUBROUTINE lbfgs_header 모형 (5) program prime 솟수 (prime number) 찾기 프로그램 implicit none integer k,kk,m do k=101,999,2 kk=k/2 do m=3,kk,2 3부터 시작하여 kk까지 변하지만, 각 단계마다 2씩 증가한다. if(k == (k/m)*m ) cycle enddo (20 of 114) 오전 10:21:34

21 write(6,*) k enddo end program prime 모형 (6) program area implicit none real*8 a,b,c,s,p,aaa read(5,*) a,b,c p=a+b+c s=p/2.0d0 aaa=sqrt(s*(s-a)*(s-b)*(s-c)) write(6,*) aaa stop end program area 모형 (7) implicit none integer n,k read(5,*) n k=2 30 if (n/k*k == n) goto 70 매우 직설적인 표현 k=k+1 if(k <= n/2) goto 30 아주 직설적으로 갈 곳을 정해주는 경우 write(6,*) n,' is a prime' stop 70 continue write(6,*) n, ' is not a prime' stop end 모형 (8) implicit none real*8 a,b,c real*8 d read(5,*) a,b,c d=b**2-4.d0*a*c if(d < 0.0d0)then write(6,*) 'NO REAL ROOTS' endif if(d == 0.0d0)then write(6,*) 'two identical roots' write(6,*) -b/(2.d0*a) endif if(d> 0.0d0)then write(6,*) 'two distinct roots' write(6,*) (-b+sqrt(d))/(2.0d0*a) (21 of 114) 오전 10:21:34

22 write(6,*) (-b-sqrt(d))/(2.0d0*a) endif stop end 모형 (9) implicit none integer n,nn real*8 a,b real(8) x,f,area,h,s integer k f(x)=((3.0d0*x-4.0d0)*x+6.0d0)*x+5.0d0 read(5,*) a,b,n h=(b-a)/float(n) nn=n-1 s=f(a)+f(b) do k=1,nn s=s+2.0d0*f(a+float(k)*h) enddo area=h*s/2.0d0 write(6,*) a,b,area stop end 모형 (10) 모형 (11) 예제 (1) 간결하면서도 강력한 포트란90의 자료처리 기법들을 맛 보시죠. 감각적이면서도 효과적인 자료처리, 마치 상대수비수들 사이로 보이는 아군에게 찔러주는 정확한 전진패스같은 시원한 맛이 나는데요. 내가 너무 오 버를 하는군요. 간결한 함수들은 컴파일러에 의해서 좀더 좋은 식 (해당 기계에 적합한 형식으로 계산이 될 것입니다. program conjugate_gradient implicit none integer iters, its, n logical :: converged real*8 :: tol, up, alpha, beta real*8, allocatable :: a(:,:), b(:), x(:), r(:), u(:), p(:), xnew(:) read(5,*) n, tol, its (22 of 114) 오전 10:21:34

23 allocate( a(n,n), b(n), x(n), r(n), u(n), p(n), xnew(n) ) 적절한 메모리 확보 open(10, file='data') read(10,*) a read(10,*) b close(10) x=1.d0 모든 컴포넌트가 1.d0로 초기화됩니다. r=b-matmul(a,x) 군더더기가 하나도 없군요. 당연히 매트릭스 곱하기 벡터죠. p=r 벡터가 복사됩니다. iters=0 do do loop 형식입니다. iters=iters+1 u=matmul(a,p) up=dot_product(r,r) 당연히 내적(inner product)이죠. alpha=up/dot_product(p,u) xnew=x+p*alpha r=r-u*alpha 여기서 알파는 숫자입니다. beta=dot_product(r,r)/up p=r+p*beta beta는 숫자입니다. converged=(maxval(abs(xnew-x))/maxval(abs(x)) < tol) 컴포넌트들 다 조사하면서 가장 큰 컴포넌트를 찾습 니다. x=xnew 당연히 벡터의 복사이지요. if( converged.or. iters == its) exit do 루프의 탈출 신호입니다. enddo write(6,*) iters write(6,*) x deallocate( a,b,x,r,u,p,xnew ) end program conjugate_gradient 이제는 필요없어요. 토사구팽 real*8 avec(10),bvec(10), avalue=sum(avec*bvec) 는 결국, avec(1)*bvec(1)+avec(2)*bvec(2)+...와 같다. 또는, avalue=dot_product(avec,bvec) 처럼 표현할 수 있다. complex 양일 경우 sum(conjg(avec)*bvec) 와 같은 응용도 가능하다. 만약, logical avec(3),bvec(3)인 경우, (avec(1).and. bvec(1)).or. (avec(2).and. bvec(2)).or...를 의미 한다. 한 가지 주의할 점: matmul(amat,bmat)는 amat*bmat와 같지 않다. 행렬 곱합기가 정의될 때, 곱한 행렬 의 ij 원소는 sum(amat(i,:)*bmat(:,j)) 와 같다. 이와 같이 어떠한 일들을 할 때 우리는 procedure를 만든다. 하지 만, 많은 경우 이미 전문가들에 의해서 잘 개발된 procedure들이 있다. 이를 이용하는 것은 매우 효율적 (일반 으로 상당히 빠른 계산을 보장한다.), 경제적일 뿐만 아니라 안정적이다. 포트란90은 113가지 instrinsic procedure가 있다. 이들을 적극 이용하는 것이 좋겠다. 뿐만 아니라, 전문가들이 개발해둔 것들도 많이 있다. BLAS (Basic Linear Algebra Subroutines) : vector, matrix-vector, matrix-matrix 계산들. LaPACK (linear algebra package) : 이들은 아주 최적화된 것이다. ATLAS, FFTW, NAG, IMSL : 참조: Python을 위 하여 개발된 모듈 Numeric python은 ATLAS를 사용한다. 예제 (2) hermitian matrix diagonalization implicit none integer nm,n,matz,ierr real*8, allocatable :: ar(:,:),ai(:,:),w(:),zr(:,:),zi(:,:),fv1(:),fv2(:),fm1(:,:) real*8 tmp integer i,j,k (23 of 114) 오전 10:21:34

24 matz=0 matz=1 nm=3 n=nm open(1,file='input',form='formatted') read(1,*) n close(1) write(6,*) n,' n' allocate(ar(nm,n),ai(nm,n),w(n),zr(nm,n),zi(nm,n),fv1(n),fv2(n),fm1(2,n)) do i=1,n do j=1,i 소위 lower-triangular-form으로 저장하고 있습니다 [ set up the hermitian matrix : lower triangular form ar(i,j)=dfloat(i+j)**2 ai(i,j)=1.d0/dfloat(i-j) ] if(i == j) ai(i,j)=0.0d0 hermitian 이기 때문에 반드시 만족해야 함 enddo enddo do i=1,n do j=1,i ar(j,i)=ar(i,j) ai(j,i)=-ai(i,j) enddo enddo call ch(nm,n,ar,ai,w,matz,zr,zi,fv1,fv2,fm1,ierr) write(6,*) 'eigenvalues' do i=1,n write(6,*) w(i) enddo write(6,*) 'eigenvectors' do j=1,n write(6,*) (dcmplx(zr(k,j),zi(k,j)),k=1,n) enddo do i=1,n do j=1,n tmp=0.0d0 do k=1,n tmp=tmp+dcmplx(zr(k,i),-zi(k,i))*dcmplx(zr(k,j),zi(k,j)) enddo write(6,*) i,j,tmp enddo enddo deallocate(ar,ai,w,zr,zi,fv1,fv2,fm1) stop end 더 이상 필요하지 않은 배열들의 메모리 할당을 풀어줍니다. (24 of 114) 오전 10:21:34

25 다운받은 ch.f 프로그램사용 (eispack) 입력 3 출력 3 n eigenvalues eigenvectors ( , ) ( , ) ( , ) ( , ) ( , ) ( , ) ( , E-002) ( , E-002) ( , ) E E E E E E 예제 (3) LBFGS루틴을 이용하여 함수를 국소 최적화하기입니다. local optimization 방법입니다. subroutine lbfgs_mini(object,posi,grad,nn) Written by In-Ho Lee, KRISS, March 3 (2003). USE sma, ONLY : sma_energy_force,tolatx,tocarx,l_pbc implicit none integer nn real*8 object,posi(nn),grad(nn) integer n_local,m,msave real*8, allocatable :: diag(:),w(:) real*8 eps,xtol,gtol,t1,t2,stpmin,stpmax integer iprint(2),iflag,icall,mp,lp,nwork integer jnatom,j logical diagco real*8, allocatable :: xyz(:,:),fxyz(:,:) real*8, allocatable :: xyz_car(:,:) The driver for LBFGS must always declare LB2 as EXTERNAL external lb2 (25 of 114) 오전 10:21:35

26 common /lb3/mp,lp,gtol,stpmin,stpmax n_local=nn m=5 ; msave=7 ; nwork=n_local*(2*msave+1)+2*msave allocate(diag(n_local),w(nwork)) iprint(1)= 1 ; iprint(2)= 0 We do not wish to provide the diagonal matrices Hk0, and therefore set DIAGCO to FALSE. diagco=.false. ; eps= 1.0d-4 ; xtol= 1.0d-16 ; icall=0 ; iflag=0 For a tight convergence : eps=1.d-7 jnatom=nn/3 allocate(xyz(jnatom,3),fxyz(jnatom,3)) allocate(xyz_car(jnatom,3)) if(l_pbc)then do j=1,jnatom xyz(j,1)=posi(3*(j-1)+1) xyz(j,2)=posi(3*(j-1)+2) xyz(j,3)=posi(3*(j-1)+3) enddo call tocarx(xyz,xyz_car) do j=1,jnatom posi(3*(j-1)+1)=xyz_car(j,1) posi(3*(j-1)+2)=xyz_car(j,2) posi(3*(j-1)+3)=xyz_car(j,3) enddo endif 20 continue do j=1,jnatom xyz_car(j,1)=posi(3*(j-1)+1) xyz_car(j,2)=posi(3*(j-1)+2) xyz_car(j,3)=posi(3*(j-1)+3) enddo if(l_pbc)then call tolatx(xyz_car,xyz) else xyz=xyz_car endif call sma_energy_force(xyz,fxyz,object) write(6,'(e14.7)') object do j=1,jnatom grad(3*(j-1)+1)=-fxyz(j,1) grad(3*(j-1)+2)=-fxyz(j,2) grad(3*(j-1)+3)=-fxyz(j,3) enddo lbfgs requires a gradient, not a force call lbfgs(n_local,m,posi,object,grad,diagco,diag,iprint,eps,xtol,w,iflag) if(iflag.le.0) go to 50 icall=icall + 1 We allow at most 2000 evaluations of Function and Gradient if(icall > 2000) go to 50 go to continue deallocate(diag,w) (26 of 114) 오전 10:21:35

27 if(l_pbc)then do j=1,jnatom xyz_car(j,1)=posi(3*(j-1)+1) xyz_car(j,2)=posi(3*(j-1)+2) xyz_car(j,3)=posi(3*(j-1)+3) enddo call tolatx(xyz_car,xyz) do j=1,jnatom posi(3*(j-1)+1)=xyz(j,1) posi(3*(j-1)+2)=xyz(j,2) posi(3*(j-1)+3)=xyz(j,3) enddo endif deallocate(xyz,fxyz) deallocate(xyz_car) return end subroutine lbfgs_header(energy,posi,grad,nn,itype) written by In-Ho Lee, KRISS, Jan USE sma, ONLY : sma_energy_force IMPLICIT NONE INTEGER nn CHARACTER*2 itype(1) REAL*8 energy,posi(1),grad(1) INTEGER n_local,m,msave REAL*8, allocatable :: diag(:),w(:) REAL*8 eps,xtol,gtol,t1,t2,stpmin,stpmax INTEGER iprint(2),iflag,icall,mp,lp,j,nwork LOGICAL diagco INTEGER iatom The driver for LBFGS must always declare LB2 as EXTERNAL external lb2 common /lb3/mp,lp,gtol,stpmin,stpmax n_local=nn m=5 ; msave=7 ; nwork=n_local*(2*msave+1)+2*msave allocate(diag(n_local),w(nwork)) iprint(1)= 1 ; iprint(2)= 0 We do not wish to provide the diagonal matrices Hk0, and therefore set DIAGCO to FALSE. diagco=.false. ; eps= 1.0d-3 ; xtol= 1.0d-16 ; icall=0 ; iflag=0 for a tight convergence : eps=1.d-7 diagco=.false. ; eps= 1.0d-7 ; xtol= 1.0d-16 ; icall=0 ; iflag=0 20 continue call sma_energy_force(itype,posi,grad,energy) do iatom=1,nn grad(iatom)=-grad(iatom) grad(iatom)=-grad(iatom) grad(iatom)=-grad(iatom) enddo write(6,'(e22.12,2x,a6)') energy,'energy' (27 of 114) 오전 10:21:35

28 call lbfgs(n_local,m,posi,energy,grad,diagco,diag,iprint,eps,xtol,w,iflag) if(iflag.le.0) go to 50 icall=icall + 1 We allow at most 2000 evaluations of Function and Gradient if(icall > 2000) go to 50 go to continue deallocate(diag,w) end subroutine lbfgs_header 예제 (4) dihedral angle 계산과 4개 원자위치에 대한 dihedral angle의 미분 module dihedral_angle Written by In-Ho Lee, KRISS, October 4 (2003). implicit none private save real*8 diangle,dadi(3),dadj(3),dadk(3),dadl(3) public :: mk_diangle,diangle,dadi,dadj,dadk,dadl contains subroutine mk_diangle(ri,rj,rk,rl) implicit none real*8 ri(3),rj(3),rk(3),rl(3) real*8 rij(3),rkj(3),rkl(3),xnorm_ij_kj,xnorm_kj_kl real*8 rij_kj(3),rkj_kl(3),rtmp(3),xnorm_kj,xnorm_mj,xnorm_nk real*8 rmj(3),rnk(3),xxr,xxi real*8 argu,qsign rij=ri-rj ; rkj=rk-rj ; rkl=rk-rl call cross(rij,rkj,rij_kj) call cross(rkj,rkl,rkj_kl) call cross(rij_kj,rkj_kl,rtmp) xnorm_ij_kj=sqrt(dot_product(rij_kj,rij_kj)) xnorm_kj_kl=sqrt(dot_product(rkj_kl,rkj_kl)) argu=dot_product(rij_kj,rkj_kl)/(xnorm_ij_kj*xnorm_kj_kl) ; argu=min(1.0d0,max(-1.0d0,argu)) qsign=1.0d0 ; if(dot_product(rkj,rtmp) < 0.0d0) qsign=-1.0d0 diangle=qsign*acos(argu) call cross(rij,rkj,rmj) call cross(rkj,rkl,rnk) xnorm_kj=sqrt(dot_product(rkj,rkj)) xnorm_mj=sqrt(dot_product(rmj,rmj)) xnorm_nk=sqrt(dot_product(rnk,rnk)) dadi(:)=rmj(:)*xnorm_kj/xnorm_mj**2 dadl(:)=-rnk(:)*xnorm_kj/xnorm_nk**2 xxr=dot_product(rij,rkj)/xnorm_kj**2-1.0d0 xxi=dot_product(rkl,rkj)/xnorm_kj**2 dadj(:)=xxr*dadi(:)-xxi*dadl(:) xxr=dot_product(rkl,rkj)/xnorm_kj**2-1.0d0 (28 of 114) 오전 10:21:35

29 xxi=dot_product(rij,rkj)/xnorm_kj**2 dadk(:)=xxr*dadl(:)-xxi*dadi(:) dadi=dadi*qsign dadj=dadj*qsign dadk=dadk*qsign dadl=dadl*qsign end subroutine mk_diangle subroutine cross(r1,r2,r3) implicit none real*8 r1(3),r2(3),r3(3) r3(1) = r1(2)*r2(3)-r1(3)*r2(2) r3(2) = r1(3)*r2(1)-r1(1)*r2(3) r3(3) = r1(1)*r2(2)-r1(2)*r2(1) end subroutine cross end module dihedral_angle =================================================================================== =================================================================================== 포트란90에서의 배열들과 관련된 메모리 관리: 결국 프로그램을 실행할 때 시스템의 메모리가 허용하는 한 프로그램을 다시 컴파일할 필요가 없도록 만들어 줍니다. integer m,n real*8, allocatable :: a(:,:) allocate(a(n,m)) a=0.0d0 생성된 배열 a는 당연히 다른 부 프로그램들에서 불려질 수 있다. 서브루틴, 함수의 인수로 사용될 수 있다. 특히, 여러 번 반복해서 불리어질 경우 해당하는 배열들의 메모리 할당여부가 불려지는 상황에 따라 달라질 수 있다. 즉, 처음 불려질 경우에는 메모리 할당이 되어있지 않는 상황이다. 이 경우 해당 배열의 allocation이 필요하다. 그 다음번 부터는 그 내용이 달라지기 때문에 메모리 할당부분은 필요가 없다. 이러한 상황을 대비 해서 프로그램할 필요가 있다. 또는 if(.not. allocated(a)) allocate(a(n,m)) 마찬가지로, 배열 a에 대한 메모리 할당이 되었었는지, 잘 모를 때는 그냥 아래와 같이 질문하고 행동을 취할 수 있다. deallocate(a) 또는 if(allocated(a)) deallocated(a) 아래의 프로그램에서는 기본적으로 사용되는 구문들을 나열했습니다. 한 번 눈으로만 확인하면 그냥 그 의 미를 알아차릴 수 있는 것들입니다. 물론, 정확한 의미의 파악을 위해서 책을 보는 것도 좋습니다 느낌표부터는 주석이 됩니다. PROGRAM xxx_1 IMPLICIT NONE INTEGER, PARAMETER :: idp=kind(0.d0) 파라미터는 프로그램 실행 중에 바뀌지 않는다는 의미입 (29 of 114) 오전 10:21:35

30 니다. CHARACTER*8 fnnd ; CHARACTER*10 fnnt 아래의 것과 같은 것입니다. CHARACTER(8) fnnd ; CHARACTER(10) fnnt LOGICAL ll1,ll2 INTEGER ii,jj,i,nn INTEGER :: kk,npt REAL*8 aa,pi REAL(8) bb,dx,xx REAL(idp) cc REAL(KIND=idp) dd REAL(KIND=idp) :: ee REAL(8) ff COMPLEX(idp) zz1,zz2,zz3 COMPLEX(idp) :: zz4 REAL(8), ALLOCATABLE :: xvector(:),xmatrix(:,:) REAL(8), ALLOCATABLE :: yvector(:) REAL(8), ALLOCATABLE :: ymatrix(:,:) REAL(8), ALLOCATABLE :: qvector(:) REAL(8) :: sma(3) REAL(8) smb(3) CALL DATE_AND_TIME(DATE=fnnd,TIME=fnnt) WRITE(6,*) ' ',' date ',fnnd,' time ',fnnt,' ' 파일이름이 input_file인 파일이 존재하는지 안 하는지를 판단한다. 즉, 아래의 명령이 실행되고 난 다음에 ll1이. true. 또는.false.인 지에 따라서 존재 여부가 알려지게 된다. INQUIRE(FILE='input_file',EXIST=ll1) pi=4.0d0*atan(1.0d0) WRITE(6,*) 'pi ',pi pi=acos(-1.0d0) WRITE(6,*) 'pi ',pi IF(ll1)THEN WRITE(6,*) 'file is present' ELSE WRITE(6,*) 'file is absent' ENDIF WRITE(6,*) idp,' idp' ll1=.true. ll2=.false. WRITE(6,*) ll1,' ll1' WRITE(6,*) ll2,' ll2' ii=0 jj=1 kk=-1 ll1= (1 == 2) ll2= (1 /= 2) IF(ii >= jj) ll1=.true. IF(ii == jj) ll2=.false. WRITE(6,*) ii,' ii' WRITE(6,*) jj,' jj' aa=0.0_idp note bb=0.0d0 cc=0.0d0 dd=0.0d0 zz1=0.0_idp zz2=cmplx(1.0_idp,2.0_idp) note (30 of 114) 오전 10:21:35

31 zz3=zz1+zz2 zz3=zz1-zz2 zz3=zz1*zz2 zz3=zz1/zz2 zz3=zz1**2 aa=real(zz3) bb=imag(zz3) cc=abs(zz3) zz4=conjg(zz3) WRITE(6,*) aa,' aa' WRITE(6,*) bb,' bb' WRITE(6,*) cc,' cc' WRITE(6,*) dd,' dd' WRITE(6,*) zz1,' zz1' WRITE(6,*) zz2,' zz2' 파일 단위 5번은 기본 입력 단위로 사용되고 파일 단위 6은 기본 출력 단위로 사용된다. 그 이외의 단위 1,,,,,,99는 아 래와 같이 사용되면 된다. 입/출력 모두 다 같은 방식으로 사용됨. OPEN(1,FILE='output1',FORM='FORMATTED') WRITE(1,'(i8)') ii WRITE(1,'(2i8)') ii,jj WRITE(1,'(i8,2x,i8)') ii,jj WRITE(1,'(i8,2x,i8,3x,a5)') ii,jj,'ii,jj' WRITE(1,*) aa,' aa' WRITE(1,*) bb,' bb' WRITE(1,*) cc,' cc' WRITE(1,*) dd,' dd' WRITE(1,*) zz1,' zz1' WRITE(1,*) zz2,' zz2' CLOSE(1) OPEN(2,FILE='output2',FORM='FORMATTED') WRITE(2,'(i8)') ii WRITE(2,'(2i8)') ii,jj CLOSE(2) sma=0.0d0 smb=1.0d0 ii=10 jj=10 nn=ii ALLOCATE(xvector(nn)) ; ALLOCATE(yvector(nn)) ALLOCATE(xmatrix(nn,nn)) ALLOCATE(ymatrix(nn,nn)) ALLOCATE(qvector(0:nn)) 0부터 시작하여 nn까지 index를 가지는 경우입니다. IF(ALLOCATED(xvector))THEN WRITE(6,*) 'xvector is allocated' WRITE(6,*) 'size ', SIZE(xvector) DO i=1,nn xvector(i)=float(i) ENDDO yvector=-xvector ENDIF WRITE(6,*) 'sum ', SUM(xvector) WRITE(6,*) 'sum ', SUM(yvector) yvector=sqrt(abs(xvector)) 다. 할당되어 있으면 작업들어 갑니다. 왠 '작업' 각 항목을 절대값으로 취한다음 그 항목에 0.5승을 취하라는 명령입니 IF(ALLOCATED(xmatrix))THEN WRITE(6,*) 'xmatrix is allocated' (31 of 114) 오전 10:21:35

32 WRITE(6,*) 'size ', SIZE(xmatrix) ENDIF xvector=0.0d0 yvector=0.0d0 xmatrix=0.0d0 ymatrix=transpose(xmatrix) ymatrix(:,1)=yvector(:) yvector=matmul(xmatrix,xvector) dd=dot_product(xvector,yvector) 매트릭스 와 벡터의 곱을 실행하라는 명령입니다. 내적입니다. inner product dd=maxval(xvector) ee=minval(yvector) dd=maxval(abs(xvector)) dd=maxval(abs(xvector-yvector))/(maxval(abs(xvector))+maxval(abs(yvector))+1.0d-8) aa=sum(xvector)/float(size(xvector)) bb=product(yvector) DEALLOCATE(xvector) ; DEALLOCATE(yvector) DEALLOCATE(xmatrix) DEALLOCATE(ymatrix) DEALLOCATE(qvector) 파일이 존재하면 지우라는 명령입니다. INQUIRE(FILE='del',EXIST=ll1) IF(LL1)THEN WRITE(6,*) 'del is present' OPEN(11,FILE='del') CLOSE(11,STATUS='DELETE') 즉, 파일이 존재한다는 것을 확인하고 지워 버리는 작업을 하고 있습니다. ENDIF aa=0.0d0 bb=1.0d0 npt= dx=(bb-aa)/float(npt-1) cc=0.0d0 DO i=1,npt xx=aa+float(i-1)*dx cc=cc+ff(xx) ENDDO cc=cc*dx write(6,*) cc,' cc' STOP END PROGRAM xxx_1 FUNCTION ff(x) REAL(8) ff 결과물로서 반환되는 자료형태를 표시하고 있슴. REAL(8) x ff=x*x RETURN END FUNCTION ff do i=1,n,1 read(5,*) aaa(i) enddo read(5,*) (aaa(i),i=1,n,1) (32 of 114) 오전 10:21:35

33 read(5,*) ((f(i,j),i=1,3),j=1,2) i=1,2,3이 먼저 실행되고 주어진 j=1상태에서 그 다음 j=2상태에서 i=1,2,3순서로 실행이 된다. -C 배열에서 허용된 인덱스의 범위를 벗어날 때 프로그램이 중단되도록 하는 옵션 real*8 a(3,2),yone(6) 은 3*2=6 개의 실수를 가리키는 이차원 배열을 정의한다. 처음 숫자는 열(row)를 가리키는 숫자로 두번째는 행(column)을 가리키는 숫자이다. 이를 굳이 1차원형태의 자료로 생각할 경우 포트란에서는 아래의 경우처 럼 취급한다. a(1,1) a(2,1) a(3,1) a(1,2) a(2,2) a(3,2) ij=3*(j-1)+i a(i,j)로 표시할 때 일차원 인덱싱(yone이라는 일차원 배열, 물론, 이 배열의 크기는 6이다.)은 결국 위와 같이 될 수 있습니다. do i=1,3 do j=1,2 ij=3*(j-1)+i yone(ij)=a(i,j) enddo enddo akvector(natom,3,nk)와 같은 배열을 특정한 형식의 1차원 배열로 바꾸는(사용자 정한 형식) 예를 아래에 표시 했다. do iatom=1,natom kp=3*nk*(iatom-1) do k=1,nk xone(kp+3*(k-1)+1)=akvector(iatom,1,k) xone(kp+3*(k-1)+2)=akvector(iatom,2,k) xone(kp+3*(k-1)+3)=akvector(iatom,3,k) enddo enddo 이러한 배열의 차원 변경은 실제 응용프로그램 작성에서 중요한 의미를 가진다. 일반으로 물리량은 다차원 배열일 때 이해하기가 쉽다. 하지만, 전산물리학에서는 이 양들이 단순한 변수일뿐이라면, 전산학에서 제시 하는 루틴들은 동등하게 취급하는 특정 1차원 배열일 가능성이 농후하다. 이 때 결국 위와 같은 조작이 유용 하다. 메인 함수에서 A(3,2)로 선언되어 사용되다가 서브루틴으로 전송될 때. subroutine abc1(a) implicit none real*8 a(3,1) subroutine abc2(a) implicit none real*8 a(3,2) 위의 둘다 정당한 선언이 된다. 물론, 두 번째가 더 정확한 선언이다. 하지만, 첫 번째 선언도 정당하다. 하지 만, 아래의 경우는 거의 에러를 포함한 코드일 것이다. (데이터의 구조를 파괴하고 있기 때문이다.) subroutine abc3(a) implicit none real*8 a(1,2) (33 of 114) 오전 10:21:35

34 가장 많이 사용될 수 있는 형태 (포트란77에서) subroutine abc4(a,lda) implicit none integer lda real*8 a(lda,1) written by In-Ho Lee, School of Physics, KIAS, March 17, 1999 implicit real*8 (a-h,o-z) integer seed(1),old(1) character*8 d_ch character*10 t_ch character*5 zonechars integer values(8) integer, external :: time call date_and_time(date=d_ch, time=t_ch) write(6,'(a8,3x,a10)') d_ch,t_ch call date_and_time(d_ch,t_ch,zonechars,values) write(6,'(a8,3x,a10,3x,a5,3x)') d_ch,t_ch,zonechars write(6,'(8i5)') values call system_clock(istime) iiss=time(itemp) write(6,*) 'seed(1)' read(5,*) seed(1) call random_seed call random_seed(size=k) write(6,*) k,' k' call random_seed(put=seed(1:k)) call random_seed(get=old(1:k)) write(6,*) seed(1),' seed(1)' write(6,*) old(1),' old(1)' sum=0.0 do ji=1,1 do ii=1,10000 call random_number(tmp) sum=sum+sin(tmp) enddo enddo write(6,*) sum,' sum' call system_clock(iftime,irate) iiff=time(itemp) write(6,*) (iiff-iiss),' sec by time' write(6,*) (iftime-istime)/float(irate),' sec by system_clock' stop end 아래의 서브루틴은 정수를 정해진 자리수하에서 문자로 변환시키는 일을 수행한다. 예를 들어 정수 10을 '0010'과 같이 변환시킨다. 4자리를 확보하여 문자로 변환된 경우이다. ################################################### ## COPYRIGHT (C) 1992 by Jay William Ponder ## ## All Rights Reserved ## (34 of 114) 오전 10:21:35

35 ################################################### ############################################################# ## ## ## subroutine numeral -- convert number to text string ## ## ## ############################################################# "numeral" converts an input integer number into the corresponding right- or left-justified text numeral number integer value of the number to be transformed string text string to be filled with corresponding numeral size on input, the minimal acceptable numeral length, if zero then output will be right justified, if nonzero then numeral is left-justified and padded with leading zeros as necessary; upon output, the number of non-blank characters in the numeral subroutine numeral (number,string,size) implicit none integer i,number,size,multi,pos integer length,minsize,len integer million,hunthou,tenthou integer thousand,hundred,tens,ones character*1 digit(0:9) character*(*) string logical right,negative data digit / '0','1','2','3','4','5','6','7','8','9' / set justification and size bounds for numeral string if (size.eq. 0) then right =.true. size = 1 else right =.false. end if minsize = size length = len(string) test the sign of the original number if (number.ge. 0) then negative =.false. else negative =.true. number = -number end if use modulo arithmetic to find place-holding digits million = number / multi = * million hunthou = (number-multi) / (35 of 114) 오전 10:21:35

36 multi = multi *hunthou tenthou = (number-multi) / multi = multi *tenthou thousand = (number-multi) / 1000 multi = multi *thousand hundred = (number-multi) / 100 multi = multi + 100*hundred tens = (number-multi) / 10 multi = multi + 10*tens ones = number - multi find the correct length to be used for the numeral if (million.ne. 0) then size = 7 else if (hunthou.ne. 0) then size = 6 else if (tenthou.ne. 0) then size = 5 else if (thousand.ne. 0) then size = 4 else if (hundred.ne. 0) then size = 3 else if (tens.ne. 0) then size = 2 else size = 1 end if size = min(size,length) size = max(size,minsize) convert individual digits to a string of numerals if (size.eq. 7) then string(1:1) = digit(million) string(2:2) = digit(hunthou) string(3:3) = digit(tenthou) string(4:4) = digit(thousand) string(5:5) = digit(hundred) string(6:6) = digit(tens) string(7:7) = digit(ones) else if (size.eq. 6) then string(1:1) = digit(hunthou) string(2:2) = digit(tenthou) string(3:3) = digit(thousand) string(4:4) = digit(hundred) string(5:5) = digit(tens) string(6:6) = digit(ones) else if (size.eq. 5) then string(1:1) = digit(tenthou) string(2:2) = digit(thousand) string(3:3) = digit(hundred) string(4:4) = digit(tens) string(5:5) = digit(ones) else if (size.eq. 4) then string(1:1) = digit(thousand) string(2:2) = digit(hundred) string(3:3) = digit(tens) (36 of 114) 오전 10:21:35

37 string(4:4) = digit(ones) else if (size.eq. 3) then string(1:1) = digit(hundred) string(2:2) = digit(tens) string(3:3) = digit(ones) else if (size.eq. 2) then string(1:1) = digit(tens) string(2:2) = digit(ones) else string(1:1) = digit(ones) end if right-justify if desired, and pad with blanks if (right) then do i = size, 1, -1 pos = length - size + i string(pos:pos) = string(i:i) end do do i = 1, length-size string(i:i) = ' ' end do else do i = size+1, length string(i:i) = ' ' end do end if return end =================================================================================== =================================================================================== 포트란과 C의 차이점 중, 다중 배열의 연속 저장방식이 다르다. 예를 들어, aa(100,100), 또는 aa[100][100]이 있을 경우 이들은 1차원 형식으로 메모리에 저장된다. 포트란의 경우, aa(1,1)과 aa(2,1),aa(3,1), aa(4,1),...과 같 은 형식으로 인접하여 저장된다. 항상 메모리에 저장된 데이터의 접근은 저장된 순서를 따를는 것이 유리하 다. 그 이유는 컴퓨터가 실제 계산을 할 떄는 cache를 사용하기 때문이다. 데이터 접근이 시간적으로 아주 빠 르게 되는 특별한 데이터들이 항상 있다. 즉, 캐쉬미스가 발생하지 않도록 해야한다. 서브루틴과 함수 비교 함수는 내장함수들처럼 사용된다. 내장함수는 instrinsic function로 표기된다. sin(), cos(), sqrt(), exp() 같은 것 들이 있다. 서브루틴은 반드시 call을 통하여 불리어질 수 있다. 함수는 항상 그 자료형을 가진다. (논리형, 정수형, 실수 형, 복소수형...) 서브루틴은 일반으로 여러 개의 인수들을 가질 수 있다. 물론, 인수를 한 개도 가지지 않는 서브루틴도 가능 하다. 인수를 통한 자료의 변화가 있을 수 있다. 물론, 없을 수도 있다. 그것은 서브루틴이 하는 일에 따라 달 라진다. 함수는 1개 이상, 여러 개의 인수를 계산하여 오직 하나의 결과만 함수 이름을 통해 넘겨준다. 주로 수식 연산 에 쓰인다. (37 of 114) 오전 10:21:35

38 기본 자료형 : 일반적으로 REAL(실수), INTEGER(정수), DOUBLE PRECISION(배정밀도 실수), COMPLEX (복소수), BOOLEAN(참, 거짓), CHARACTER(문자) 등의 6가지를 지원한다. 기본연산 ** : 지수 연산으로 A**2는 A 2 를 의미한다. * : 곱셈 연산 / : 나눗셈 연산 +, - : 덧셈, 뺄셈 논리 연산에는.NOT.,.AND.,.OR.가 있다 관계 연산은 크기를 판정한다..GT. : > ( 보다 더 큰가?) >.GE. : >= ( 과 크거나 같은가?) >=.LT. : < ( 보다 더 작은가?) <.LE. : <= ( 과 작거나 같은가?) <=.EQ. : = ( 과 같은가?) ==.NE. : NOT = ( 과 다른가?) /= 포트란은 6가지의 수형(type)을 가지고 있다. 포트란77 중심의 경우 아래와 같이 분류가능하다. 포트란90과 비교하면 아주 편협된 개념을 사용하고 있다. 물론, 이들만 가지고도 거의 모든일들을 프로그램할 수 있다. 문제는 프로그램의 효율적인 작성이다. 1 자료형 : 일반적으로 REAL(실수), INTEGER(정수), DOUBLE PRECISION(배정밀도 실수), COMPLEX(복 소수), BOOLEAN(참, 거짓), CHARACTER(문자) 등의 6가지를 지원한다. 2 변수(variable) : 변수는 첫 글자가 I 와 N 사이의 글자로 시작하면 정수형을 의미하며 나머지 영문자로 시 작하는 변수는 실수형이 된다. 하지만 별도로 수형을 선언하면 원래 의미가 무시된다. 변수는 영문자와 숫자 를 이용하여 6자리까지 가능하지만 첫 글자는 반드시 영문자여야 한다. 특수문자는 변수명에 사용할 수 없 다. 3 배열(array) : 배열이란 여러 개의 변수가 하나로 묶인 것을 의미하는데 배열의 첨자는 상수나 정수형 변 수 또는 상수와 정수형 변수의 조합 수식이 가능하다. 4 COMMON(공용) : 단위가 다른 프로그램에 있는 변수나 배열명 등이 같은 기억장소를 사용할 수 있는 선 언문이다. COMMON /named_common1/ A,B,C COMMON /abc_name/ AREA(5), AVRG COMMON B(10), YY, XX c_name이라는 named common block을 선언한 경우: 내용은 실수 두 개를 순서대로 배열한 것이다. (일반으로 (38 of 114) 오전 10:21:35

39 는 다양한 형태의 배열,변수들을 나열할 수 있다.) common block에 사용된 배열, 변수들은 서브루틴, 함수의 인수로 사용되서는 안된다. common block에 속하는 변수들의 순서가 중요하다. 즉, 이름은 다를 수 있다는 것 이다. 물론, 갯수는 항상 같아야 한다. 같은 장소에 순서대로 저장된 변수들로 볼 수 있다. 특정 서브루틴에서 어떤 변수를 참조할 때 저장된 순서를 보고 그 해당하는 값을 찾아오는 것이다. 예를 들어, 실수형 변수 3개 의 나열은 실수형 배열로 크기 3인 것으로 대치될 수 있다. a,b,c ==> abc(3) (위험하기도 하지만, 유용한 측면도 있다.) 시간이 지난다음에 다시 프로그램을 볼 때는 다소 위험한 것이 사 실이다. 어떤 변수는 여러 개의 named common block에 동시에 등록될 수도 있다. 실제로 많이 사용되는 형식 은 common bolck에 해당하는 변수들을 한 파일에 저장한 다음 이를 해당 common block이 사용되는 프로그램 들에서 불러들이는 방법이다. 이 때 include 'file_name.h'형식으로 불러 들인다. 물론 선언문 부분에서 불러들 이면 된다. 이렇게 하면 사용하고 있는 common block이 같은 형식이라는 것을 보장받게 된다. 파일이 하나이 기 때문에 항상 같은 형식이 필요한 루틴에 복사가 되는 것이다. 아주 강력한 방법이다. subroutine subx (...) integer nmax parameter (nmax=20) integer n real A(nmax, nmax) common /matrix/ A, n, nmax nmax의 값은 main program 뿐만 아니라 위의 block이 사용되는 모든 곳에서 그 변수들의 종류와 크기가 정확 하게 일치하여야 한다. 또한 matrix의 크기는 compile하기 전에 알려져야 하기에 nmax는 parameter 문장에서 정의하여야만 한다. 위에서 표시된 부분을 하나의 파일에 저장한다. 예를 들어 matrix.h라는 파일에. 그 다음 에는 이를 필요로 하는 모든 곳에서 불러들인다. include 'matrix.h'형태로 말이다. subroutine subx (...) include 'matrix.h' matrix size를 변경할 때는 위의 파일을 수정하고 다시 컴파일하면 될 것이다. 이런 방식을 사용하지 않을 경 우, 위의 block이 사용되는 모든 곳을 찾아다니면서 똑같은 사이즈를 맞추어 줘야 한다. 대단히 불편하며 오 류로 이어질 위험성이 짙다. include를 사용한다 하더라도 문제는 있을 수 있다. 변수들의 이름이 우연히 같 을 수 있기 때문이다. 한 쪽에서 잘 정의된 변수명이 다른 쪽 프로그램에서 하필이면 다른 뜻으로 사용될 수 있기 때문이다. 꼭 체크를 해야 한다. (이것이 common block 사용의 약점이다. 이를 강제로 하는 방법은 없을 까?) program main real*8 alpha, beta real*8 a,b,aa,bb integer iii common /c_name/ alpha, beta common /dd_name/ a,b,iii,aa,bb stop end subroutine sub1 () real*8 alpha, beta common /c_name/ alpha, beta (39 of 114) 오전 10:21:35

40 return end subroutine sub2 () integer jjj real*8 a,b,aabb(2) real*8 alpha1,beta1 common /c_name/ alpha1, beta1 common /dd_name/ a,b,jjj,aabb(2) return end 5 EQUIVALENCE(동등) : 하나의 단위 프로그램 내에서 변수나 배열 등이 같은 기억장소를 공유할 수 있게 하는 선언문이다. 사용하는 변수의 이름은 달라도 같은곳에 저장됨을 의미한다. common block의 경우와 마 찬가지이다. 포트란 90을 배우는데 필요없는 것을 논할 필요가 있을까? 하지만, 옛날 소스를 분석할 때를 위 해서 알아둘 필요가 있다. 6 DATA문 : 초기값을 부여할 수 있다. 아래의 예를 보면, 배열 A와 변수 PI, X가 나타나 있다. PI에는 초기 값 3.14가 배정되며 X에는 1.0이 배정된다. 배열 A는 3개 모두에 1.0이 배정된다. 예1) DATA문 DATA A,B,C/10.5, 20.1,30.2/ 예2) DATA문 DIMENSION M(25) DATA M/25*4/ 예3) 암시적 DO를 사용한 DATA문 DIMENSION A(5,6), B(10) DATA (A(I,6), I=1,5)/5*3.5) DATA (B(J),J=1,10)/5*2.1, 3*5.0, 7.0, 8.0/ DATA KG, LG, MG /10, 20, 30/ DIMENSION A(20) DATA A /20*8/ KG=10, LG=20, MG=30, 배열 A의 모든 원소에 8로 초기와 됨 배열을 선언하는 다른 방법이 있다. real*8 A, x dimension x(-50:50) dimension A(0:10,-2:20) 은 real*8 A(0:10,-2:20), x(-50:50) 과 같다. =================================================================================== (40 of 114) 오전 10:21:35

41 =================================================================================== 유닉스 명령어 make에 대해서 설명드리겠습니다. 프로그램 개발과 사용에 아주 유용한 명령어입니다. Introduction to Make makefile이라는 파일명이 있습니다. 연관되어 사용되는 유닉스 명령어는 make입니다. 한 디렉토리에 모여있는 (한곳이 아니라도 상관은 없습니다.) 소스(원시코드)들을 한꺼번에 컴파일하는 것 을 포함하여 이들을 쉽게 관리하는 것입니다. 예를 들어 하나의 실행화일 만들 때 또는 개발 중에 있는 코드들 중에서 하나가 수정되어 바뀌었을 때 (에러 를 잡는다. 또는 기능을 추가한다. 등등..) 유사한 일반적인 경우에 자동으로 연관된 프로그램들만 컴파일 해서 실행화일을 다시 만들 수 있는 명령어가 make입니다. 기본적으로 makefile이 존재하는 디렉토리에서 make라는 유닉스 명령어를 치면 makefile안에 있는 명령들이 실행됩니다. 어떠한 명령이 수행되는가는 makefile을 만드는 방식에 따라서 달라집니다. 본란에서는 아주 기 초적인 것만 다루려고 합니다. 가장 간단한 실행은 아마도 "순서대로 컴파일을 수행하여 실행화일을 만듭니다." 일 것입니다. 하나의 루틴이상이 소스 수준에서 고쳐졌다면, 사용자가 make라는 명령어를 사용하면 자동으로 "변화가 있 는 소스들을 알아서 컴파일을 하여 새로운 실행화일을 만들어 줍니다. 신기한 기능이지요. 날짜를 체크해서 변화가 있는지를 알아 내는 것입니다. (실제 유닉스에서는) 예를 들면, 아래의 예에서 알 수 있듯이 amd.x라는 실행화일을 만든다고 합시다. 모든 object 프로그램들이 필 요할 것입니다. 그래서 변수 이름을 AMD라고 두고 있습니다.당연히 AMD=.선언하는 것이죠. "\"는 다 음 줄로 확장할 때 상용하는 문자입니다. 명령어가 한 줄을 넘어간다고 가정하자. 이때 그냥 줄줄이 적는다 면 읽기도 힘들고, 작성하는 사람도 조금 찜찜하다. 이때 '\' 문자를 이용해서 여러 라인으로 나타낼 수 있다. 변수가 사용될 때는 $()처럼 사용됩니다. 아래의 문법(2줄에 대해서 설명에서는 두 번째 줄에서 시작할 때 tab이 사용됨을 꼭 숙지하시기 바람 니다. amd.x라는 실행화일은 $(AMD)로 선언된 것들로부터 만들어진다. 만드는 방법은 두 번째줄에 있다. 이 때 tab을 이용하여 다음줄을 시작한다. 구체적으로 컴파일러를 이용하 여 컴파일하는 것을 보여줌. 이 경우 오브젝트들을 연결하는 경우임. 메인을 만들기 때문에. target : dependency command amd.x: $(AMD) $(segl) -o amd.x $(AMD) $(LINK) $(LIB) 아래의 경우를 설명합니다. 이 경우 오브젝트 energy_force.o가 어떻게 연결되어 있는가를 알려 줍니다. 소위 dependence라고 하지요. 예를 들어서 tight_binding_lapack_c.f90 라는 프로그램을 조금 고쳤습니다. 그러면 이 함수(모듈, 서브루틴, 함수)를 실제로 이용하고 있는 energy_force.f90 루틴도 새로이 컴파일되어야 합니다. (41 of 114) 오전 10:21:35

42 이런 것들을 dependence라고 하지요. 이를 정확히 가려쳐 주면 편리합니다. 당연히, energy_force.o 오브젝트는 energy_force.f90라는 것으로 만들고, 여기에 use를 사용하여 모듈 tight_binding_lapack_c.f90를 이용하고 있는 경우는 아래와 같이 그 dependence를 알려줄 필요가 있습니다. energy_force.o:energy_force.f90 tight_binding_lapack_c.o $(cmpl) $(Ocft) energy_force.f90 "탭 (tab)"이 들어간 곳의 위치를 확인해야 한다. 두 줄단위로 복사해서 사용하면 tab을 잊어 버리지 않을 수 있게 된다 cmpl=pgf90 segl=pgf90 Ocft= -c -fast cft= -c #LINK=-L/usr/lib/xmm -L/usr/lib/xmm/atlas LINK= -L/usr/lib/sse2 LIBS = -llapack -lf77blas -latlas -lf2c AMD= amd.o predictor_corrector_nose.o energy_force.o tight_binding_lapack_c.o \ mini_driver.o minimization.o lbfgs.o bond_boost.o erf.o maxwell_boltzmann.o \ check_transition.o rmarin.o angle_boost.o dihedral.o dihedral_boost.o amd.x: $(AMD) $(segl) -o amd.x $(AMD) $(LINK) $(LIBS) amd.o:amd.f90 predictor_corrector_nose.o bond_boost.o angle_boost.o dihedral_boost.o $(cmpl) $(Ocft) amd.f90 energy_force.o:energy_force.f90 tight_binding_lapack_c.o $(cmpl) $(Ocft) energy_force.f90 tight_binding_lapack_c.o:tight_binding_lapack_c.f90 $(cmpl) $(Ocft) tight_binding_lapack_c.f90 predictor_corrector_nose.o:predictor_corrector_nose.f90 $(cmpl) $(Ocft) predictor_corrector_nose.f90 mini_driver.o:mini_driver.f90 $(cmpl) $(Ocft) mini_driver.f90 minimization.o:minimization.f90 $(cmpl) $(Ocft) minimization.f90 bond_boost.o:bond_boost.f90 $(cmpl) $(Ocft) bond_boost.f90 angle_boost.o:angle_boost.f90 $(cmpl) $(Ocft) angle_boost.f90 dihedral_boost.o:dihedral_boost.f90 dihedral.o $(cmpl) $(Ocft) dihedral_boost.f90 dihedral.o:dihedral.f90 erf.o:erf.f $(cmpl) $(Ocft) erf.f clean: rm -f *.x *.o *.mod *.M core* touch: touch *.f90 *.i makefile ; chmod 600 *.f90 *.i makefile ; ls -l *.f90 *.i makefile rmo: rm -f *.o *.mod *.M core* lsl: ls -l *.f90 makefile *.i (42 of 114) 오전 10:21:35

43 소소들이 있는 곳에 두고 사용자의 유닉스 환경에 따라서 조금 수정을 하십시요. 특히 tab이 있는 곳을 주의 하십시오. 소스들의 이름을 지금 가지고 있는 것들로 바꾸고 마찬가지로 오브젝트들의 이름도 마찬가지입니 다. dependece를 확실히 체크하여 주십시오. pgf90을 ifc로 바꾸고 다 바꾼 다음 가능한 명령어는 첫 번째는 오브젝트를 모두 다 지우는 것입니다. (make clean) 두 번째는 실행 화일을 만드는 것입니다. (make) 또는 소스를 약간 수정한 후 바로 make하면 됩니다. 그렇게 되면 최종적으 로 바뀐 것만 다시 컴파일해서 실행화일을 만듭니다. 프로그램 개발 단계에서 아주 유용하겠죠. 개발 후에도 물론 유용하구요. touch는 유닉스 명령어입니다. 파일들의 생성시간을 이 명령어에 의해서 재조정하게 합니 다. 재조정되는 시간은 '터치되는 시간입니다. 추후 어떠한 파일들이 변화가 되었는지를 체크할 수 있게 해줍 니다. 예를 들어, 유닉스 명령어 ls에서 추가 옵션을 사용하면 됩니다. ls -ltra 는 변화된 파일들의 시간 순서상 으로 파일들을 보여줍니다. ls -Fs는 디렉토리 표시와 파일 사이즈 크기 표시 make clean make edit : main.o kbd.o command.o display.o \ insert.o search.o files.o utils.o cc -o edit main.o kbd.o command.o display.o \ insert.o search.o files.o utils.o main.o : main.c defs.h cc -c main.c kbd.o : kbd.c defs.h command.h cc -c kbd.c command.o : command.c defs.h command.h cc -c command.c display.o : display.c defs.h buffer.h cc -c display.c insert.o : insert.c defs.h buffer.h cc -c insert.c search.o : search.c defs.h buffer.h cc -c search.c files.o : files.c defs.h buffer.h command.h cc -c files.c utils.o : utils.c defs.h cc -c utils.c clean : rm edit main.o kbd.o command.o display.o \ insert.o search.o files.o utils.o 이 `edit'라고 불리는 실행 파일을 만들도록 makefile을 쓰기 위해서는 다음과 같이 입력한다: (43 of 114) 오전 10:21:35

44 make 이 makefile을 사용해서 실행 파일과 모든 오브젝트 파일들을 삭제하려면 다음과 같이 입력한다: make clean 이 예제 makefile에서 target들은 `edit'라는 실행 파일과 `main.o'와 `kbd.o'라는 오브젝트 파일들을 포함한다. dependencies는 `main.c'와 `defs.h'와 같은 파일들이다. 실제로 각 `.o' 파일은 target이면서 동시에 dependency이 다. command는 `cc -c main.c'와 `cc -c kbd.c'를 포함한다 target이 파일일 때, 어떤 dependencies 들 중 하나가 변경되었을 때 재 컴파일되거나 재 링크되어야 한다. 게다 가 스스로 자동적으로 생성된 어떤 dependencies는 맨 먼저 업데이트되어야 한다. 이 예제에서 `edit'는 8개의 오브젝트 파일들에 의존한다; `main.o'는 소스 파일 `main.c'와 헤더 파일 `defs.h'에 의존한다. objects = main.o kbd.o command.o display.o \ insert.o search.o files.o utils.o edit : $(objects) cc -o edit $(objects) main.o : main.c defs.h cc -c main.c kbd.o : kbd.c defs.h command.h cc -c kbd.c command.o : command.c defs.h command.h cc -c command.c display.o : display.c defs.h buffer.h cc -c display.c insert.o : insert.c defs.h buffer.h cc -c insert.c search.o : search.c defs.h buffer.h cc -c search.c files.o : files.c defs.h buffer.h command.h cc -c files.c utils.o : utils.c defs.h cc -c utils.c clean : rm edit $(objects) objects = main.o kbd.o command.o display.o \ insert.o search.o files.o utils.o edit : $(objects) (44 of 114) 오전 10:21:35

45 cc -o edit $(objects) main.o : defs.h kbd.o : defs.h command.h command.o : defs.h command.h display.o : defs.h buffer.h insert.o : defs.h buffer.h search.o : defs.h buffer.h files.o : defs.h buffer.h command.h utils.o : defs.h.phony : clean clean : -rm edit $(objects) objects = main.o kbd.o command.o display.o \ insert.o search.o files.o utils.o edit : $(objects) cc -o edit $(objects) $(objects) : defs.h kbd.o command.o files.o : command.h display.o insert.o search.o files.o : buffer.h =================================================================================== =================================================================================== netcdf (network Common Data Format) 컴퓨터 기종에 관계없는 하나의 과학적 자료저장 형식을 말합니다. 1. self-describing 2. portable 3. direct acess (저장된 자료들 중에서 필요한 부분에 직접접근한다. <--> sequential access) 4. appendable 5. sharable 포트란77을 위한 netcdf 인터페이스 라이버러리들이 있었고, 많이 사용되었다. 최근 포트란90을 직접지원하는 인터페이스 라이브러리들은 상당히 간결한 형태로 netcdf를 지원한다. Hierarchical Data Format (HDF) (45 of 114) 오전 10:21:35

46 예제 (2) ####################################################################### Code converted using TO_F90 by Alan Miller Date: Time: 17:25:32 MODULE GA_commons IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) include 'params.f' INTEGER, PARAMETER :: indmax=200, nchrmax=30, nparmax=2 COMMON /ga1/ npopsiz, nowrite INTEGER, SAVE :: npopsiz, nowrite COMMON /ga2/ nparam, nchrome INTEGER, SAVE :: nparam, nchrome COMMON /ga3/ parent, iparent REAL (dp), SAVE :: parent(nparmax,indmax) INTEGER, SAVE :: iparent(nchrmax,indmax) COMMON /ga4/ fitness REAL (dp), SAVE :: fitness(indmax) COMMON /ga5/ g0, g1, ig2 REAL (dp), SAVE :: g0(nparmax), g1(nparmax) INTEGER, SAVE :: ig2(nparmax) COMMON /ga6/ parmax, parmin, pardel, nposibl REAL (dp), SAVE :: parmax(nparmax), parmin(nparmax), pardel(nparmax) INTEGER, SAVE :: nposibl(nparmax) COMMON /ga7/ child, ichild REAL (dp), SAVE :: child(nparmax,indmax) INTEGER, SAVE :: ichild(nchrmax,indmax) COMMON /ga8/ nichflg INTEGER, SAVE :: nichflg(nparmax) COMMON /inputga/ pcross, pmutate, pcreep, maxgen, idum, irestrt, & itourny, ielite, icreep, iunifrm, iniche, iskip, iend, nchild, & microga, kountmx REAL (dp), SAVE :: pcross, pmutate, pcreep INTEGER, SAVE :: maxgen, idum, irestrt, itourny, ielite, icreep, & iunifrm, iniche, iskip, iend, nchild, microga, kountmx END MODULE GA_commons (46 of 114) 오전 10:21:35

47 예제 (3) random 함수 MODULE random A module for random number generation from the following distributions: Distribution Function/subroutine name Normal (Gaussian) random_normal Gamma random_gamma Chi-squared random_chisq Exponential random_exponential Weibull random_weibull Beta random_beta t random_t Multivariate normal random_mvnorm Generalized inverse Gaussian random_inv_gauss Poisson random_poisson Binomial random_binomial1 * random_binomial2 * Negative binomial random_neg_binomial von Mises random_von_mises Cauchy random_cauchy Generate a random ordering of the integers 1.. N random_order Initialize (seed) the uniform random number generator for ANY compiler seed_random_number ** Two functions are provided for the binomial distribution. If the parameter values remain constant, it is recommended that the first function is used (random_binomial1). If one or both of the parameters change, use the second function (random_binomial2). The compilers own random number generator, SUBROUTINE RANDOM_NUMBER(r), is used to provide a source of uniformly distributed random numbers. N.B. At this stage, only one random number is generated at each call to one of the functions above. The module uses the following functions which are included here: bin_prob to calculate a single binomial probability lngamma to calculate the logarithm to base e of the gamma function Some of the code is adapted from Dagpunar's book: Dagpunar, J. 'Principles of random variate generation' Clarendon Press, Oxford, ISBN In most of Dagpunar's routines, there is a test to see whether the value of one or two floating-point parameters has changed since the last call. These tests have been replaced by using a logical variable FIRST. This should be set to.true. on the first call using new values of the parameters, and.false. if the parameter values are the same as for the previous call. (47 of 114) 오전 10:21:35

48 Version 1.10, 1 August 1996 Changes from version The random_order, random_poisson & random_binomial routines have been replaced with more efficient routines. 2. A routine, seed_random_number, has been added to seed the uniform random number generator. This requires input of the required number of seeds for the particular compiler from a specified I/O unit such as a keyboard. Author: Alan Miller CSIRO Division of Mathematics & Statistics Private Bag 10, Clayton South MDC Clayton 3169, Victoria, Australia Phone: (+61) Fax: (+61) mel.dms.csiro.au WWW-page: REAL, PRIVATE :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0, & vsmall = TINY(1.0), vlarge = HUGE(1.0) PRIVATE :: integral CONTAINS REAL FUNCTION random_normal() Adapted from the following Fortran 77 code ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM. THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOL. 18, NO. 4, DECEMBER, 1992, PP The function random_normal() returns a normally distributed pseudo-random number with zero mean and unit variance. The algorithm uses the ratio of uniforms method of A.J. Kinderman and J.F. Monahan augmented with quadratic bounding curves. IMPLICIT NONE Local variables REAL :: s = , t = , a = , b = , & r1 = , r2 = , u, v, x, y, q Generate P = (u,v) uniform in rectangle enclosing acceptance region DO CALL RANDOM_NUMBER(u) CALL RANDOM_NUMBER(v) v = * (v - half) Evaluate the quadratic form x = u - s y = ABS(v) - t q = x**2 + y*(a*y - b*x) Accept P if inside inner ellipse (48 of 114) 오전 10:21:35

49 IF (q < r1) EXIT Reject P if outside outer ellipse IF (q > r2) CYCLE Reject P if outside acceptance region IF (v**2 < -4.0*LOG(u)*u**2) EXIT Return ratio of P's coordinates as the normal deviate random_normal = v/u RETURN END FUNCTION random_normal REAL FUNCTION random_gamma(s, first) Adapted from Fortran 77 code from the book: Dagpunar, J. 'Principles of random variate generation' Clarendon Press, Oxford, ISBN FUNCTION GENERATES A RANDOM GAMMA VARIATE. CALLS EITHER random_gamma1 (S > 1.0) OR random_exponential (S = 1.0) OR random_gamma2 (S < 1.0). S = SHAPE PARAMETER OF DISTRIBUTION (0 < REAL). IMPLICIT NONE REAL, INTENT(IN) :: s LOGICAL, INTENT(IN) :: first IF (s <= zero) THEN PRINT *, 'SHAPE PARAMETER VALUE MUST BE POSITIVE' STOP IF (s > one) THEN random_gamma = random_gamma1(s, first) ELSE IF (s < one) THEN random_gamma = random_gamma2(s, first) ELSE random_gamma = random_exponential() RETURN END FUNCTION random_gamma REAL FUNCTION random_gamma1(s, first) Adapted from Fortran 77 code from the book: Dagpunar, J. 'Principles of random variate generation' Clarendon Press, Oxford, ISBN FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO GAMMA**(S-1)*EXP(-GAMMA), BASED UPON BEST'S T DISTRIBUTION METHOD S = SHAPE PARAMETER OF DISTRIBUTION (1.0 < REAL) (49 of 114) 오전 10:21:35

50 IMPLICIT NONE REAL, INTENT(IN) :: s LOGICAL, INTENT(IN) :: first Local variables REAL :: b, h, d, r, g, f, x, sixty4 = 64.0, three = 3.0, pt75 = 0.75 SAVE b, h IF (s <= one) THEN PRINT *, 'IMPERMISSIBLE SHAPE PARAMETER VALUE' STOP IF (first) THEN Initialization, if necessary b = s - one h = SQRT(three*s - pt75) DO CALL RANDOM_NUMBER(r) g = r - r*r IF (g <= zero) CYCLE f = (r - half)*h/sqrt(g) x = b + f IF (x <= zero) CYCLE CALL RANDOM_NUMBER(r) d = sixty4*g*(r*g)**2 IF (d <= zero) EXIT IF (d*x < x - two*f*f) EXIT IF (LOG(d) < two*(b*log(x/b) - f)) EXIT random_gamma1 = x RETURN END FUNCTION random_gamma1 REAL FUNCTION random_gamma2(s, first) Adapted from Fortran 77 code from the book: Dagpunar, J. 'Principles of random variate generation' Clarendon Press, Oxford, ISBN FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO GAMMA2**(S-1) * EXP(-GAMMA2), USING A SWITCHING METHOD. S = SHAPE PARAMETER OF DISTRIBUTION (REAL < 1.0) IMPLICIT NONE REAL, INTENT(IN) :: s LOGICAL, INTENT(IN) :: first Local variables REAL :: a, p, c, uf, vr, d, r, x, w (50 of 114) 오전 10:21:35

51 SAVE a, p, c, uf, vr, d IF (s <= zero.or. s >= one) THEN PRINT *, 'SHAPE PARAMETER VALUE OUTSIDE PERMITTED RANGE' STOP IF (first) THEN Initialization, if necessary a = one - s p = a/(a + s*exp(-a)) IF (s < vsmall) THEN PRINT *, 'SHAPE PARAMETER VALUE TOO SMALL' STOP c = one/s uf = p*(vsmall/a)**s vr = one - vsmall d = a*log(a) DO CALL RANDOM_NUMBER(r) IF (r >= vr) THEN CYCLE ELSE IF (r > p) THEN x = a - LOG((one - r)/(one - p)) w = a*log(x)-d ELSE IF (r > uf) THEN x = a*(r/p)**c w = x ELSE random_gamma2 = zero RETURN CALL RANDOM_NUMBER(r) IF (one-r <= w.and. r > zero) THEN IF (r*(w + one) >= one) CYCLE IF (-LOG(r) <= w) CYCLE EXIT random_gamma2 = x RETURN END FUNCTION random_gamma2 REAL FUNCTION random_chisq(ndf, first) Generates a random variate from the chi-squared distribution with ndf degrees of freedom IMPLICIT NONE INTEGER, INTENT(IN) :: ndf LOGICAL, INTENT(IN) :: first random_chisq = two * random_gamma(half*float(ndf), first) RETURN (51 of 114) 오전 10:21:35

52 END FUNCTION random_chisq REAL FUNCTION random_exponential( ) Adapted from Fortran 77 code from the book: Dagpunar, J. 'Principles of random variate generation' Clarendon Press, Oxford, ISBN FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL TO EXP(-random_exponential), USING INVERSION. IMPLICIT NONE Local variable REAL :: r DO CALL RANDOM_NUMBER(r) IF (r > zero) EXIT random_exponential = -LOG(r) RETURN END FUNCTION random_exponential REAL FUNCTION random_weibull(a) Generates a random variate from the Weibull distribution with probability density: a a-1 -x f(x) = a.x e IMPLICIT NONE REAL, INTENT(IN) :: a For speed, there is no checking that a is not zero or very small. random_weibull = random_exponential() ** (one/a) RETURN END FUNCTION random_weibull REAL FUNCTION random_beta(aa, bb, first) Adapted from Fortran 77 code from the book: Dagpunar, J. 'Principles of random variate generation' Clarendon Press, Oxford, ISBN FUNCTION GENERATES A RANDOM VARIATE IN [0,1] FROM A BETA DISTRIBUTION WITH DENSITY (52 of 114) 오전 10:21:35

53 PROPORTIONAL TO BETA**(AA-1) * (1-BETA)**(BB-1). USING CHENG'S LOG LOGISTIC METHOD. AA = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL) BB = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL) IMPLICIT NONE REAL, INTENT(IN) :: aa, bb LOGICAL, INTENT(IN) :: first Local variables REAL, PARAMETER :: aln4 = REAL :: a, b, d, f, h, t, c, g, r, s, x, y, z LOGICAL :: swap SAVE d, f, h, t, c, swap IF (aa <= zero.or. bb <= zero) THEN PRINT *, 'IMPERMISSIBLE SHAPE PARAMETER VALUE(S)' STOP IF (first) THEN Initialization, if necessary a = aa b = bb swap = b > a IF (swap) THEN g = b b = a a = g d = a/b f = a+b IF (b > one) THEN h = SQRT((two*a*b - f)/(f - two)) t = one ELSE h = b t = one/(one + (a/(vlarge*b))**b) c = a+h DO CALL RANDOM_NUMBER(r) CALL RANDOM_NUMBER(x) s = r*r*x IF (r < vsmall.or. s <= zero) CYCLE IF (r < t) THEN x = LOG(r/(one - r))/h y = d*exp(x) z = c*x + f*log((one + d)/(one + y)) - aln4 IF (s - one > z) THEN IF (s - s*z > one) CYCLE IF (LOG(s) > z) CYCLE random_beta = y/(one + y) ELSE IF (4.0*s > (one + one/d)**f) CYCLE random_beta = one EXIT (53 of 114) 오전 10:21:35

54 IF (swap) random_beta = one - random_beta RETURN END FUNCTION random_beta REAL FUNCTION random_t(m) Adapted from Fortran 77 code from the book: Dagpunar, J. 'Principles of random variate generation' Clarendon Press, Oxford, ISBN FUNCTION GENERATES A RANDOM VARIATE FROM A T DISTRIBUTION USING KINDERMAN AND MONAHAN'S RATIO METHOD. M = DEGREES OF FREEDOM OF DISTRIBUTION (1 <= 1NTEGER) IMPLICIT NONE INTEGER, INTENT(IN) :: m Local variables REAL :: s, c, a, f, g, r, x, v, three = 3.0, four = 4.0, quart = 0.25, & five = 5.0, sixteen = 16.0 INTEGER :: mm = 0 SAVE s, c, a, f, g IF (m < 1) THEN PRINT *, 'IMPERMISSIBLE DEGREES OF FREEDOM' STOP IF (m.ne. mm) THEN Initialization, if necessary s = m c = -quart*(s + one) a = four/(one + one/s)**c f = sixteen/a IF (m > 1) THEN g = s - one g = ((s + one)/g)**c*sqrt((s+s)/g) ELSE g = one mm = m DO CALL RANDOM_NUMBER(r) IF (r <= zero) CYCLE CALL RANDOM_NUMBER(v) x = (two*v - one)*g/r v = x*x IF (v > five - a*r) THEN IF (m >= 1.AND. r*(v + three) > f) CYCLE IF (r > (one + v/s)**c) CYCLE EXIT (54 of 114) 오전 10:21:35

55 random_t = x RETURN END FUNCTION random_t SUBROUTINE random_mvnorm(n, h, d, f, first, x, ier) Adapted from Fortran 77 code from the book: Dagpunar, J. 'Principles of random variate generation' Clarendon Press, Oxford, ISBN N.B. An extra argument, ier, has been added to Dagpunar's routine SUBROUTINE GENERATES AN N VARIATE RANDOM NORMAL VECTOR USING A CHOLESKY DECOMPOSITION. ARGUMENTS: N = NUMBER OF VARIATES IN VECTOR (INPUT,INTEGER >= 1) H(J) = J'TH ELEMENT OF VECTOR OF MEANS (INPUT,REAL) X(J) = J'TH ELEMENT OF DELIVERED VECTOR (OUTPUT,REAL) D(J*(J-1)/2+I) = (I,J)'TH ELEMENT OF VARIANCE MATRIX (J> = I) (INPUT,REAL) F((J-1)*(2*N-J)/2+I) = (I,J)'TH ELEMENT OF LOWER TRIANGULAR DECOMPOSITION OF VARIANCE MATRIX (J <= I) (OUTPUT,REAL) FIRST =.TRUE. IF THIS IS THE FIRST CALL OF THE ROUTINE OR IF THE DISTRIBUTION HAS CHANGED SINCE THE LAST CALL OF THE ROUTINE. OTHERWISE SET TO.FALSE. (INPUT,LOGICAL) ier = 1 if the input covariance matrix is not +ve definite = 0 otherwise IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL, INTENT(IN) :: h(n), d(n*(n+1)/2) REAL, INTENT(OUT) :: f(n*(n+1)/2), x(n) LOGICAL, INTENT(IN) :: first INTEGER, INTENT(OUT) :: ier Local variables INTEGER :: n2, j, i, m REAL :: y, v SAVE n2 IF (n < 1) THEN PRINT *, 'SIZE OF VECTOR IS NON POSITIVE' STOP ier = 0 IF (first) THEN Initialization, if necessary n2 = 2*n (55 of 114) 오전 10:21:35

56 IF (d(1) < zero) THEN ier = 1 RETURN f(1) = SQRT(d(1)) y = one/f(1) DO j = 2,n f(j) = d(1+j*(j-1)/2) * y DO i = 2,n v = d(i*(i-1)/2+i) DO m = 1,i-1 v = v - f((m-1)*(n2-m)/2+i)**2 IF (v < zero) THEN ier = 1 RETURN v = SQRT(v) y = one/v f((i-1)*(n2-i)/2+i) = v DO j = i+1,n v = d(j*(j-1)/2+i) DO m = 1,i-1 v = v - f((m-1)*(n2-m)/2+i)*f((m-1)*(n2-m)/2 + j) m = 1,i-1 f((i-1)*(n2-i)/2 + j) = v*y j = i+1,n i = 2,n x(1:n) = h(1:n) DO j = 1,n y = random_normal() DO i = j,n x(i) = x(i) + f((j-1)*(n2-j)/2 + i) * y i = j,n j = 1,n END SUBROUTINE random_mvnorm REAL FUNCTION random_inv_gauss(h, b, first) Adapted from Fortran 77 code from the book: Dagpunar, J. 'Principles of random variate generation' Clarendon Press, Oxford, ISBN FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY] FROM A REPARAMETERISED GENERALISED INVERSE GAUSSIAN (GIG) DISTRIBUTION WITH DENSITY PROPORTIONAL TO GIG**(H-1) * EXP(-0.5*B*(GIG+1/GIG)) USING A RATIO METHOD. H = PARAMETER OF DISTRIBUTION (0 <= REAL) B = PARAMETER OF DISTRIBUTION (0 < REAL) IMPLICIT NONE REAL, INTENT(IN) :: h, b (56 of 114) 오전 10:21:35

57 LOGICAL, INTENT(IN) :: first Local variables REAL :: a, c, d, e, ym, xm, r, w, r1, r2, x, quart = 0.25 SAVE a, c, d, e IF (h < zero.or. b <= zero) THEN PRINT *, 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES' STOP IF (first) THEN Initialization, if necessary IF (h > quart*b*sqrt(vlarge)) THEN PRINT *, 'THE RATIO H:B IS TOO SMALL' STOP e = b*b d = h + one ym = (-d + SQRT(d*d + e))/b IF (ym < vsmall) THEN PRINT *, 'THE VALUE OF B IS TOO SMALL' STOP d = h - one xm = (d + SQRT(d*d + e))/b d = half*d e = -quart*b r = xm + one/xm w = xm*ym a = w**(-half*h) * SQRT(xm/ym) * EXP(-e*(r - ym - one/ym)) IF (a < vsmall) THEN PRINT *, 'THE VALUE OF H IS TOO LARGE' STOP c = -d*log(xm) - e*r DO CALL RANDOM_NUMBER(r1) IF (r1 <= zero) CYCLE CALL RANDOM_NUMBER(r2) x = a*r2/r1 IF (x <= zero) CYCLE IF (LOG(r1) < d*log(x) + e*(x + one/x) + c) EXIT random_inv_gauss = x RETURN END FUNCTION random_inv_gauss INTEGER FUNCTION random_poisson(mu, first) ********************************************************************** Translated to Fortran 90 by Alan Miller from: RANLIB Library of Fortran Routines for Random Number Generation Compiled and Written by: Barry W. Brown (57 of 114) 오전 10:21:35

58 James Lovato Department of Biomathematics, Box 237 The University of Texas, M.D. Anderson Cancer Center 1515 Holcombe Boulevard Houston, TX This work was supported by grant CA from the National Cancer Institute. GENerate POIsson random deviate Function Generates a single random deviate from a Poisson distribution with mean mu. Arguments mu --> The mean of the Poisson distribution from which a random deviate is to be generated. REAL mu Method For details see: Ahrens, J.H. and Dieter, U. Computer Generation of Poisson Deviates From Modified Normal Distributions. ACM Trans. Math. Software, 8, 2 (June 1982), IMPLICIT NONE TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL SEPARATION OF CASES A AND B.. Scalar Arguments.. REAL, INTENT(IN) :: mu LOGICAL, INTENT(IN) :: first.... Local Scalars.. REAL a0, a1, a2, a3, a4, a5, a6, a7, b1, b2, c, c0, c1, c2, c3, d, del, & difmuk, e, fk, fx, fy, g, omega, p, p0, px, py, q, s, t, u, v, x, xx INTEGER j, k, kflag, l, m LOGICAL :: full_init.... Local Arrays.. REAL fact(10), pp(35).... Intrinsic Functions.. INTRINSIC ABS, LOG, EXP, FLOAT, IFIX, MAX, MIN, SIGN, SQRT.... Data statements.. DATA a0, a1, a2, a3, a4, a5, a6, a7/-.5, , , , & , , , / DATA fact/1., 1., 2., 6., 24., 120., 720., 5040., , / (58 of 114) 오전 10:21:35

59 SAVE pp, s, d, l, m, p, q, p0, full_init.... Executable Statements.. IF (mu > 10.0) THEN C A S E A. (RECALCULATION OF S, D, L IF MU HAS CHANGED) IF (first) THEN s = SQRT(mu) d = 6.0*mu*mu THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU ) IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10. l = IFIX(mu ) full_init =.false. STEP N. NORMAL SAMPLE - random_normal() FOR STANDARD NORMAL DEVIATE g = mu + s*random_normal() IF (g > 0.0) THEN random_poisson = IFIX(g) STEP I. IMMEDIATE ACCEPTANCE IF random_poisson IS LARGE ENOUGH IF (random_poisson>=l) RETURN STEP S. SQUEEZE ACCEPTANCE - SAMPLE U fk = FLOAT(random_Poisson) difmuk = mu - fk CALL RANDOM_NUMBER(u) IF (d*u >= difmuk*difmuk*difmuk) RETURN STEP P. PREPARATIONS FOR STEPS Q AND H. (RECALCULATIONS OF PARAMETERS IF NECESSARY) =(2*PI)**(-.5) E-1=1./ =1./7. THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK. C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION. IF (.NOT. full_init) THEN omega = /s b1 = E-1/mu b2 =.3*b1*b1 c3 = *b1*b2 c2 = b2-15.*c3 c1 = b1-6.*b *c3 c0 = 1. - b1 + 3.*b2-15.*c3 c =.1069/mu full_init =.true. IF (g<0.0) GO TO 50 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN) (59 of 114) 오전 10:21:35

60 kflag = 0 GO TO 70 STEP Q. QUOTIENT ACCEPTANCE (RARE CASE) 40 IF (fy-u*fy <= py*exp(px-fx)) RETURN STEP E. EXPONENTIAL SAMPLE - random_exponential() FOR STANDARD EXPONENTIAL DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT' (IF T <= THEN PK < FK FOR ALL MU >= 10.) 50 e = random_exponential() CALL RANDOM_NUMBER(u) u = u + u - one t = SIGN(e, u) IF (t <= (-.6744)) GO TO 50 random_poisson = IFIX(mu+s*t) fk = FLOAT(random_Poisson) difmuk = mu - fk 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN) kflag = 1 GO TO 70 STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION) 60 IF (c*abs(u) > py*exp(px+e) - fy*exp(fx+e)) GO TO 50 RETURN STEP F. 'SUBROUTINE' F. CALCULATION OF PX, PY, FX, FY. CASE random_poisson < 10 USES FACTORIALS FROM TABLE FACT 70 IF (random_poisson>=10) GO TO 80 px = -mu py = mu**random_poisson/fact(random_poisson+1) GO TO 110 CASE random_poisson >= 10 USES POLYNOMIAL APPROXIMATION A0-A7 FOR ACCURACY WHEN ADVISABLE E-1=1./ =(2*PI)**(-.5) 80 del = E-1/fk del = del - 4.8*del*del*del v = difmuk/fk IF (ABS(v)>0.25) THEN px = fk*log(one + v) - difmuk - del ELSE px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) - del py = /SQRT(fk) 110 x = (half - difmuk)/s xx = x*x fx = -half*xx fy = omega* (((c3*xx + c2)*xx + c1)*xx + c0) IF (kflag) 40, 40, C A S E B. mu < 10 (60 of 114) 오전 10:21:35

61 START NEW TABLE AND CALCULATE P0 IF NECESSARY ELSE IF (first) THEN m = MAX(1, IFIX(mu)) l = 0 p = EXP(-mu) q = p p0 = p STEP U. UNIFORM SAMPLE FOR INVERSION METHOD DO CALL RANDOM_NUMBER(u) random_poisson = 0 IF (u <= p0) RETURN STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE PP-TABLE OF CUMULATIVE POISSON PROBABILITIES (0.458=PP(9) FOR MU=10) IF (l == 0) GO TO 150 j = 1 IF (u > 0.458) j = MIN(l, m) DO k = j, l IF (u <= pp(k)) GO TO 180 IF (l == 35) CYCLE STEP C. CREATION OF NEW POISSON PROBABILITIES P AND THEIR CUMULATIVES Q=PP(K) 150 l = l + 1 DO k = l, 35 p = p*mu/float(k) q = q + p pp(k) = q IF (u <= q) GO TO 170 l = l = k 180 random_poisson = k RETURN END FUNCTION random_poisson INTEGER FUNCTION random_binomial1(n, p, first) FUNCTION GENERATES A RANDOM BINOMIAL VARIATE USING C.D.Kemp's method. This algorithm is suitable when many random variates are required with the SAME parameter values for n & p. P = BERNOULLI SUCCESS PROBABILITY (0 <= REAL <= 1) N = NUMBER OF BERNOULLI TRIALS (1 <= INTEGER) FIRST =.TRUE. for the first call using the current parameter (61 of 114) 오전 10:21:35

62 values =.FALSE. if the values of (n,p) are unchanged from last call Reference: Kemp, C.D. (1986). `A modal method for generating binomial variables', Commun. Statist. - Theor. Meth. 15(3), IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL, INTENT(IN) :: p LOGICAL, INTENT(IN) :: first Local variables INTEGER :: r0, ru, rd REAL :: odds_ratio, u, p_r, pd, pu, zero = 0.0, one = 1.0 SAVE odds_ratio, r0, p_r IF (first) THEN r0 = (n+1)*p p_r = bin_prob(n, p, r0) odds_ratio = p / (one - p) CALL RANDOM_NUMBER(u) u = u - p_r IF (u < zero) THEN random_binomial1 = r0 RETURN pu = p_r ru = r0 pd = p_r rd = r0 DO rd = rd - 1 IF (rd >= 0) THEN pd = pd * (rd+1) / (odds_ratio * FLOAT(n-rd)) u = u - pd IF (u < zero) THEN random_binomial1 = rd RETURN ru = ru + 1 IF (ru <= n) THEN pu = pu * (n-ru+1) * odds_ratio / FLOAT(ru) u = u - pu IF (u < zero) THEN random_binomial1 = ru RETURN This point should not be reached, but just in case: random_binomial1 = r0 (62 of 114) 오전 10:21:35

63 RETURN END FUNCTION random_binomial1 REAL FUNCTION bin_prob(n, p, r) Calculate a binomial probability IMPLICIT NONE INTEGER, INTENT(IN) :: n, r REAL, INTENT(IN) :: p Local variable REAL :: one = 1.0 bin_prob = EXP( lngamma(dble(n+1)) - lngamma(dble(r+1)) - lngamma(dble(n-r+1)) & + r*log(p) + (n-r)*log(one - p) ) RETURN END FUNCTION bin_prob DOUBLE PRECISION FUNCTION lngamma(x) Logarithm to base e of the gamma function. Accurate to about 1.e-14. Programmer: Alan Miller Latest revision of Fortran 77 version - 28 February 1988 IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: x Local variables DOUBLE PRECISION :: a1 = D-02, a2 = D-03, & a3 = D-04, a4 = D-04, & temp, arg, product, lnrt2pi = D-1, & pi = D0 LOGICAL :: reflect lngamma is not defined if x = 0 or a negative integer. IF (x > 0.d0) GO TO 10 IF (x /= INT(x)) GO TO 10 lngamma = 0.d0 RETURN If x < 0, use the reflection formula: gamma(x) * gamma(1-x) = pi * cosec(pi.x) 10 reflect = (x < 0.d0) IF (reflect) THEN arg = 1.d0 - x ELSE arg = x (63 of 114) 오전 10:21:35

64 Increase the argument, if necessary, to make it > 10. product = 1.d0 20 IF (arg <= 10.d0) THEN product = product * arg arg = arg + 1.d0 GO TO 20 Use a polynomial approximation to Stirling's formula. N.B. The real Stirling's formula is used here, not the simpler, but less accurate formula given by De Moivre in a letter to Stirling, which is the one usually quoted. arg = arg - 0.5D0 temp = 1.d0/arg**2 lngamma = lnrt2pi + arg * (LOG(arg) - 1.d0 + & (((a4*temp + a3)*temp + a2)*temp + a1)*temp) - LOG(product) IF (reflect) THEN temp = SIN(pi * x) lngamma = LOG(pi/temp) - lngamma RETURN END FUNCTION lngamma INTEGER FUNCTION random_binomial2(n, pp, first) ********************************************************************** Translated to Fortran 90 by Alan Miller from: RANLIB Library of Fortran Routines for Random Number Generation Compiled and Written by: Barry W. Brown James Lovato Department of Biomathematics, Box 237 The University of Texas, M.D. Anderson Cancer Center 1515 Holcombe Boulevard Houston, TX This work was supported by grant CA from the National Cancer Institute. GENerate BINomial random deviate Function Generates a single random deviate from a binomial distribution whose number of trials is N and whose probability of an event in each trial is P. Arguments N --> The number of trials in the binomial distribution from which a random deviate is to be generated. INTEGER N (64 of 114) 오전 10:21:35

65 P --> The probability of an event in each trial of the binomial distribution from which a random deviate is to be generated. REAL P FIRST --> Set FIRST =.TRUE. for the first call to perform initialization the set FIRST =.FALSE. for further calls using the same pair of parameter values (N, P). LOGICAL FIRST random_binomial2 <-- A random deviate yielding the number of events from N independent trials, each of which has a probability of event P. INTEGER random_binomial Method This is algorithm BTPE from: Kachitvichyanukul, V. and Schmeiser, B. W. Binomial Random Variate Generation. Communications of the ACM, 31, 2 (February, 1988) 216. ********************************************************************** IMPLICIT NONE *****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY.... Scalar Arguments.. REAL, INTENT(IN) :: pp INTEGER, INTENT(IN) :: n LOGICAL, INTENT(IN) :: first.... Local Scalars.. REAL :: al, alv, amaxp, c, f, f1, f2, ffm, fm, g, p, p1, p2, p3, p4, q, & qn, r, u, v, w, w2, x, x1, x2, xl, xll, xlr, xm, xnp, xnpq, xr, & ynorm, z, z2 REAL :: zero = 0.0, half = 0.5, one = 1.0 INTEGER :: i, ix, ix1, k, m, mp.... Intrinsic Functions.. INTRINSIC ABS, LOG, MIN, INT, SQRT SAVE p, q, xnp, ffm, m, fm, xnpq, p1, xm, xl, xr, c, al, xll, xlr, & p2, p3, p4, qn, r, g.... Executable Statements.. *****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE IF (first) THEN p = MIN(pp, one-pp) q = one - p xnp = n * p (65 of 114) 오전 10:21:35

66 IF (xnp > 30.) THEN IF (first) THEN ffm = xnp + p m = ffm fm = m xnpq = xnp * q p1 = INT(2.195*SQRT(xnpq) - 4.6*q) + half xm = fm + half xl = xm - p1 xr = xm + p1 c = / ( fm) al = (ffm-xl) / (ffm - xl*p) xll = al * (one + half*al) al = (xr - ffm) / (xr*q) xlr = al * (one + half*al) p2 = p1 * (one + c + c) p3 = p2 + c / xll p4 = p3 + c / xlr *****GENERATE VARIATE, Binomial mean at least CALL RANDOM_NUMBER(u) u = u * p4 CALL RANDOM_NUMBER(v) TRIANGULAR REGION IF (u <= p1) THEN ix = xm - p1 * v + u GO TO 110 PARALLELOGRAM REGION IF (u <= p2) THEN x = xl + (u-p1) / c v = v * c + one - ABS(xm-x) / p1 IF (v > one.or. v <= zero) GO TO 20 ix = x ELSE LEFT TAIL IF (u <= p3) THEN ix = xl + LOG(v) / xll IF (ix < 0) GO TO 20 v = v * (u-p2) * xll ELSE RIGHT TAIL ix = xr - LOG(v) / xlr IF (ix > n) GO TO 20 v = v * (u-p3) * xlr *****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST (66 of 114) 오전 10:21:35

67 k = ABS(ix-m) IF (k <= 20.OR. k >= xnpq/2-1) THEN EXPLICIT EVALUATION f = one r = p / q g = (n+1) * r IF (m < ix) THEN mp = m + 1 DO i = mp, ix f = f * (g/i-r) ELSE IF (m > ix) THEN ix1 = ix + 1 DO i = ix1, m f = f / (g/i-r) IF (v > f) THEN GO TO 20 ELSE GO TO 110 SQUEEZING USING UPPER AND LOWER BOUNDS ON LOG(F(X)) amaxp = (k/xnpq) * ((k*(k/ ) )/xnpq + half) ynorm = -k * k / (2.*xnpq) alv = LOG(v) IF (alv<ynorm - amaxp) GO TO 110 IF (alv>ynorm + amaxp) GO TO 20 STIRLING'S (actually de Moivre's) FORMULA TO MACHINE ACCURACY FOR THE FINAL ACCEPTANCE/REJECTION TEST x1 = ix + 1 f1 = fm + one z = n fm w = n - ix + one z2 = z * z x2 = x1 * x1 f2 = f1 * f1 w2 = w * w IF (alv - (xm*log(f1/x1) + (n-m+half)*log(z/w) + (ix-m)*log(w*p/(x1*q)) + & ( (462.-(132.-( /f2)/f2)/f2)/f2)/f1/ & ( (462.-(132.-( /z2)/z2)/z2)/z2)/z/ & ( (462.-(132.-( /x2)/x2)/x2)/x2)/x1/ & ( (462.-(132.-( /w2)/w2)/w2)/w2)/w/ ) > zero) THEN GO TO 20 ELSE GO TO 110 ELSE INVERSE CDF LOGIC FOR MEAN LESS THAN 30 IF (first) THEN qn = q ** n (67 of 114) 오전 10:21:35

68 r = p / q g = r * (n+1) 90 ix = 0 f = qn CALL RANDOM_NUMBER(u) 100 IF (u >= f) THEN IF (ix > 110) GO TO 90 u = u - f ix = ix + 1 f = f * (g/ix - r) GO TO IF (pp > half) ix = n - ix random_binomial2 = ix RETURN END FUNCTION random_binomial2 INTEGER FUNCTION random_neg_binomial(sk, p) Adapted from Fortran 77 code from the book: Dagpunar, J. 'Principles of random variate generation' Clarendon Press, Oxford, ISBN FUNCTION GENERATES A RANDOM NEGATIVE BINOMIAL VARIATE USING UNSTORED INVERSION AND/OR THE REPRODUCTIVE PROPERTY. SK = NUMBER OF FAILURES REQUIRED (Dagpunar's words) = the `power' parameter of the negative binomial (0 < REAL) P = BERNOULLI SUCCESS PROBABILITY (0 < REAL < 1) THE PARAMETER H IS SET SO THAT UNSTORED INVERSION ONLY IS USED WHEN P <= H, OTHERWISE A COMBINATION OF UNSTORED INVERSION AND THE REPRODUCTIVE PROPERTY IS USED. IMPLICIT NONE REAL, INTENT(IN) :: sk, p Local variables THE PARAMETER ULN = -LOG(MACHINE'S SMALLEST REAL NUMBER). REAL, PARAMETER :: h = 0.7 REAL :: q, x, st, uln, v, r, s, y, g INTEGER :: k, i, n IF (sk <= zero.or. p <= zero.or. p >= one) THEN PRINT *, 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES' STOP q = one - p x = zero st = sk (68 of 114) 오전 10:21:35

69 IF (p > h) THEN v = one/log(p) k = st DO i = 1,k DO CALL RANDOM_NUMBER(r) IF (r > zero) EXIT n = v*log(r) x = x + n st = st - k s = zero uln = -LOG(vsmall) IF (st > -uln/log(q)) THEN PRINT *, ' P IS TOO LARGE FOR THIS VALUE OF SK' STOP y = q**st g = st CALL RANDOM_NUMBER(r) DO IF (y > r) EXIT r = r - y s = s + one y = y*p*g/s g = g + one random_neg_binomial = x + s + half RETURN END FUNCTION random_neg_binomial REAL FUNCTION random_von_mises(k, first, ier) Algorithm VMD from: Dagpunar, J.S. (1990) `Sampling from the von Mises distribution via a comparison of random numbers', J. of Appl. Statist., 17, Fortran 90 code by Alan Miller CSIRO Division of Mathematics & Statistics Arguments: k (real) parameter of the von Mises distribution. first (logical) set to.true. the first time that the function is called, or the first time with a new value for k. When first =.TRUE., the function sets up starting values and may be very much slower. ier (integer) error indicator = 0 if k is acceptable = -1 if k < 0. = +1 if k > = 10. IMPLICIT NONE REAL, INTENT(IN) :: k LOGICAL, INTENT(IN) :: first INTEGER, INTENT(OUT) :: ier (69 of 114) 오전 10:21:35

70 Local variables INTEGER :: nk, j, n REAL, PARAMETER :: pi = REAL :: p(20), theta(0:20), sump, r, th, lambda, rlast DOUBLE PRECISION :: dk SAVE nk, p, theta IF (first) THEN Initialization, if necessary ier = 0 IF (k < zero) THEN ier = -1 RETURN nk = k + k + one IF (nk > 20) THEN ier = +1 RETURN dk = k theta(0) = zero IF (k > half) THEN Set up array p of probabilities. sump = zero DO j = 1, nk IF (j < nk) THEN theta(j) = ACOS(one - j/k) ELSE theta(nk) = pi Numerical integration of e^[k.cos(x)] from theta(j-1) to theta(j) CALL integral(theta(j-1), theta(j), p(j), dk) sump = sump + p(j) j = 1, nk p(1:nk) = p(1:nk) / sump ELSE p(1) = one theta(1) = pi if k > 0.5 if first CALL RANDOM_NUMBER(r) DO j = 1, nk r = r - p(j) IF (r < zero) EXIT r = -r/p(j) DO th = theta(j-1) + r*(theta(j) - theta(j-1)) lambda = k - j + one - k*cos(th) n = 1 rlast = lambda (70 of 114) 오전 10:21:35

71 DO CALL RANDOM_NUMBER(r) IF (r > rlast) EXIT n = n + 1 rlast = r IF (n.ne. 2*(n/2)) EXIT is n even? CALL RANDOM_NUMBER(r) random_von_mises = SIGN(th, (r - rlast)/(one - rlast) - half) RETURN END FUNCTION random_von_mises SUBROUTINE integral(a, b, result, dk) Gaussian integration of exp(k.cosx) from a to b. IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: dk REAL, INTENT(IN) :: a, b REAL, INTENT(OUT) :: result Local variables DOUBLE PRECISION :: xmid, range, x1, x2, & x(3) = (/ D0, D0, D0/), & w(3) = (/ D0, D0, D0/) INTEGER :: i xmid = (a + b)/2.d0 range = (b - a)/2.d0 result = 0.d0 DO i = 1, 3 x1 = xmid + x(i)*range x2 = xmid - x(i)*range result = result + w(i)*(exp(dk*cos(x1)) + EXP(dk*COS(x2))) result = result * range RETURN END SUBROUTINE integral REAL FUNCTION random_cauchy() Generate a random deviate from the standard Cauchy distribution IMPLICIT NONE Local variables REAL :: v(2) DO (71 of 114) 오전 10:21:35

72 CALL RANDOM_NUMBER(v) v = two*(v - half) IF (ABS(v(2)) < vsmall) CYCLE Test for zero IF (v(1)**2 + v(2)**2 < one) EXIT random_cauchy = v(1) / v(2) END FUNCTION random_cauchy SUBROUTINE random_order(order, n) Generate a random ordering of the integers 1... n. IMPLICIT NONE INTEGER, INTENT(IN) :: n INTEGER, INTENT(OUT) :: order(n) Local variables INTEGER :: i, j, k REAL :: wk DO i = 1, n order(i) = i Starting at the end, swap the current last indicator with one randomly chosen from those preceeding it. DO i = n, 2, -1 CALL RANDOM_NUMBER(wk) j = 1 + i * wk IF (j < i) THEN k = order(i) order(i) = order(j) order(j) = k RETURN END SUBROUTINE random_order SUBROUTINE seed_random_number(iounit) IMPLICIT NONE INTEGER, INTENT(IN) :: iounit Local variables INTEGER :: k INTEGER, ALLOCATABLE :: seed(:) CALL RANDOM_SEED(SIZE=k) ALLOCATE( seed(k) ) WRITE(*, '(1x, a, i2, a)')'enter ', k, ' integers for random no. seeds: ' (72 of 114) 오전 10:21:35

73 READ(*, *) seed WRITE(iounit, '(1x, a, (7i10))') 'Random no. seeds: ', seed CALL RANDOM_SEED(PUT=seed) DEALLOCATE( seed ) END SUBROUTINE seed_random_number END MODULE random 예제 (4) global optimization MODULE global_minimum This is translated from the code of: Dr. Tibor Csendes Dept. of Applied Informatics, Jozsef Attila University H-6701 Szeged, Pf. 652, Hungary Phone: Fax: [email protected] URL: IMPLICIT NONE INTEGER, PARAMETER, PRIVATE :: dp = SELECTED_REAL_KIND(12, 60) CONTAINS ROUTINE NAME - GLOBAL Code converted using TO_F90 by Alan Miller Date: Time: 13:16: LATEST REVISION - OKTOBER 23, 1986 Latest revision of Fortran 90 version - 30 January 2000 PURPOSE - GLOBAL MINIMUM OF FUNCTION OF N VARIABLES USING A LOCAL SEARCH METHOD USAGE - CALL GLOBAL (AMIN, AMAX, NPARM, M, N100, NG0, IPR, NSIG, X0, NC, F0) ARGUMENTS AMIN - VECTOR OF LENGTH NPARM CONTAINING THE LOWER BOUNDS OF THE PARAMETERS, SO X(I) IS SEARCHED IN THE INTERVAL (AMIN(I), AMAX(I)). (INPUT) AMAX - VECTOR OF LENGTH NPARM CONTAINING THE UPPER BOUNDS OF THE PARAMETERS. (INPUT) NPARM - NUMBER OF PARAMETERS, 1 <= NPARM <= 20. (INPUT) (73 of 114) 오전 10:21:35

74 M - NUMBER OF RESIDUAL FUNCTIONS, WHEN THE OBJECTIVE FUNCTION IS OF THE FORM F1**2 + F2** FM**2, 1 <=M <= 100. (INPUT) N.B. M is NOT used for this purpose It is passed to the user's routine FUNCT as a parameter. N100 - NUMBER OF SAMPLE POINTS TO BE DRAWN UNIFORMLY IN ONE CYCLE, 20 <= N100 <= THE SUGGESTED VALUE IS 100*NPARM. (INPUT) NG0 - NUMBER OF BEST POINTS SELECTED FROM THE ACTUAL SAMPLE, 1 <= NG0 <= 20. THE SUGGESTED VALUE IS TWICE THE EXPECTED NUMBER OF LOCAL MINIMA. (INPUT) IPR - FORTRAN DATA SET REFERENCE NUMBER WHERE THE PRINTED OUTPUT BE SENT. (INPUT) NSIG - CONVERGENCE CRITERION, THE ACCURACY REQUIRED IN THE PARAMETER ESTIMATES. THIS CONVERGENCE CRITERION IS SATISFIED IF ON TWO SUCCESSIVE ITERATIONS THE PARAMETER ESTIMATES AGREE, COMPONENT BY COMPONENT, TO NSIG DIGITS. THE SUGGESTED VALUE IS 6. (INPUT) X0 - OUTPUT NPARM BY 20 MATRIX CONTAINING NC (UP TO 20) LOCAL MINIMIZERS FOUND. NC - NUMBER OF DIFFERENT LOCAL MINIMIZERS FOUND. (OUTPUT) F0 - OUTPUT VECTOR OF NC (UP TO 20) OBJECTIVE FUNCTION VALUES, F0(I) BELONGS TO THE PARAMETERS X0(1,I), X0(2,I),..., X0(NPARM,I). REQUIRED ROUTINES - URDMN, FUN, LOCAL SUBROUTINE global(amin, amax, nparm, m, n100, ng0, ipr, nsig, x0, nc, f0) REAL (dp), INTENT(IN) :: amin(:) REAL (dp), INTENT(IN) :: amax(:) INTEGER, INTENT(IN) :: nparm INTEGER, INTENT(IN) :: m INTEGER, INTENT(IN OUT) :: n100 INTEGER, INTENT(IN OUT) :: ng0 INTEGER, INTENT(IN) :: ipr INTEGER, INTENT(IN) :: nsig REAL (dp), INTENT(OUT) :: x0(:,:) INTEGER, INTENT(OUT) :: nc REAL (dp), INTENT(OUT) :: f0(:) REAL (dp) :: x(nparm,100), x1(nparm,20), xcl(nparm,100), r(nparm), w(nparm) INTEGER :: ic(100), ic1(20) INTEGER :: i, i1, icc, icj, ig, ii, iii, im, in1, inum, inum1, inum2, it, & iv, j, jj, l1, maxfn, n, n0, n1, ncp, nfe, nfe1, ng, ng10, nm, & nn100, ns REAL (dp) :: f(100), f1(20), fcl(100), y(nparm), mmin(nparm), mmax(nparm) REAL (dp) :: a, alfa, b, b1, bb, fc, ff, fm, relcon REAL (dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp, two = 2.0_dp, ten = 10._dp DO i=1,nparm mmin(i) = amin(i) mmax(i) = amax(i) IF (mmin(i) == mmax(i)) GO TO 460 b1 = one/nparm IF (ng0 < 1) ng0 = 1 (74 of 114) 오전 10:21:35

75 IF (ng0 > 20) ng0 = 20 IF (n100 < 20) n100 = 20 IF (n100 > 10000) n100 = IF (n100 >= 100) GO TO 10 nn100 = n100 n = 1 GO TO nn100 = 100 n = n100/100 n100 = n* ng10 = 100 DO i=1,ng10 f(i) = HUGE(one) ic(i) = 0 DO i=1,nparm mmax(i) = (mmax(i) - mmin(i))/two mmin(i) = mmin(i) + mmax(i) alfa = 0.01_dp nfe = 0 ng = 0 ns = 0 nc = 0 ncp = 1 n0 = 0 n1 = 0 im = 1 ig = 0 fm = HUGE(one) maxfn = 500*nparm relcon = ten**(-nsig) SAMPLING 20 n0 = n0 + n100 nm = n0-1 ng = ng + ng0 ns = ns + 1 IF (ns*ng0 > 100) GO TO 465 b = (one - alfa**(one/real(nm)))**b1 bb = 0.1*b DO i1=1,n DO j=1,nn100 CALL RANDOM_NUMBER( r ) DO i=1,nparm y(i) = two*r(i) - one CALL fun(y, fc, nparm, m, mmin, mmax) IF (fc >= fm) CYCLE f(im) = fc DO i=1,nparm x(i,im) = y(i) IF (im <= ng.and. ic(im) > 0) ig = ig - 1 ic(im) = 0 im = 1 fm = f(1) DO i=2,ng10 IF (f(i) < fm) CYCLE (75 of 114) 오전 10:21:35

76 im = i fm = f(i) nfe = nfe + n100 WRITE(ipr,901) n100 WRITE(*,901) n FORMAT(/' ', i5, ' FUNCTION EVALUATIONS USED FOR SAMPLING') SORTING inum = ng10-1 DO i=1,inum im = i fm = f(i) inum1 = i + 1 DO j=inum1,ng10 IF (f(j) >= fm) CYCLE im = j fm = f(j) IF (im <= i) CYCLE a = fm DO j=1,nparm y(j) = x(j,im) IF (i > ng.or. im <= ng) GO TO 55 IF (ic(ng) == 0.AND. ic(im) > 0) ig = ig+1 IF (ic(ng) > 0.AND. ic(im) == 0) ig = ig-1 55 icc = ic(im) inum1 = im-i DO j=1,inum1 inum2 = im-j f(inum2+1) = f(inum2) ic(inum2+1) = ic(inum2) DO jj=1,nparm x(jj,inum2+1) = x(jj,inum2) f(i) = a DO j=1,nparm x(j,i) = y(j) ic(i) = icc IF (nc <= 0) GO TO 200 CLUSTERING TO X* DO iii=1,nc i = 1 in1 = i fcl(i) = f0(iii) DO j=1,nparm xcl(j,i) = x0(j,iii) DO j=1,ng IF (ic(j) /= iii) CYCLE in1 = in1+1 DO ii=1,nparm xcl(ii,in1) = x(ii,j) (76 of 114) 오전 10:21:35

77 95 DO j=1,ng IF (ic(j) /= 0) CYCLE IF (fcl(i) >= f(j)) CYCLE DO l1=1,nparm w(l1) = ABS(xcl(l1,i)-x(l1,j)) a = zero DO l1=1,nparm IF (w(l1) > a) a = w(l1) IF (a >= b) CYCLE WRITE(ipr,902) iii WRITE(*,902) iii 902 FORMAT(' SAMPLE POINT ADDED TO THE CLUSTER NO. ', i2) DO ii=1,nparm w(ii) = x(ii,j)*mmax(ii) + mmin(ii) WRITE(ipr,903) f(j), w(1:nparm) WRITE (*,903) f(j), w(1:nparm) 903 FORMAT(' ', g14.8/ (' ', 5(g14.8, ' '))) ig = ig+1 IF (ig >= ng) GO TO 395 in1 = in1+1 fcl(in1) = f(j) DO ii=1,nparm xcl(ii,in1) = x(ii,j) ic(j) = iii i = i+1 IF (i <= in1) GO TO 95 IF (n1 <= 0) GO TO 200 CLUSTERING TO X1 DO iii=1,n1 i = 1 in1 = i fcl(i) = f1(iii) DO j=1,nparm xcl(j,i) = x1(j,iii) 155 DO j=1,ng IF (ic(j) /= 0) CYCLE IF (fcl(i) >= f(j)) CYCLE DO l1=1,nparm w(l1) = ABS(xcl(l1,i)-x(l1,j)) a = zero DO l1=1,nparm IF (w(l1) > a) a=w(l1) IF (a >= b) CYCLE WRITE(ipr,902) ic1(iii) WRITE(*,902) ic1(iii) DO ii=1,nparm w(ii) = x(ii,j)*mmax(ii) + mmin(ii) WRITE(ipr,903) f(j), w(1:nparm) WRITE(*,903) f(j), w(1:nparm) (77 of 114) 오전 10:21:35

78 ig = ig+1 IF (ig >= ng) GO TO 395 in1 = in1+1 fcl(in1) = f(j) DO ii=1,nparm xcl(ii,in1) = x(ii,j) ic(j) = ic1(iii) i = i+1 IF (i <= in1) GO TO 155 LOCAL SEARCH 200 it = 0 DO i1=1,ng IF (ic(i1) /= 0) CYCLE DO i=1,nparm y(i) = x(i,i1) ff = f(i1) CALL local(m, nparm, relcon, maxfn, y, ff, nfe1, mmin, mmax) IF (nc <= 0) GO TO 290 DO iv=1,nc DO l1=1,nparm w(l1) = ABS(x0(l1,iv) - y(l1)) a = zero DO l1=1,nparm IF (w(l1) > a) a = w(l1) IF (a < bb) GO TO 255 GO TO 290 NEW SEED-POINT 255 n1 = n1 + 1 WRITE(ipr,905) iv, nfe1 WRITE(*,905) iv,nfe1 905 FORMAT(' NEW SEED POINT ADDED TO THE CLUSTER NO. ', i2, ', NFEV=', i5) DO ii=1,nparm w(ii) = x(ii,i1)*mmax(ii) + mmin(ii) WRITE(ipr,903) ff, w(1:nparm) WRITE(*,903) ff, w(1:nparm) IF (ff >= f0(iv)) GO TO 280 WRITE(ipr,906) iv, f0(iv), ff WRITE(*,906) iv,f0(iv),ff 906 FORMAT(' *** IMPROVEMENT ON THE LOCAL MINIMUM NO. ', & i2, ':', g14.8, ' FOR ', g14.8) w(1:nparm) = y(1:nparm)*mmax(1:nparm) + mmin(1:nparm) WRITE(ipr,903) ff, w(1:nparm) WRITE(*,903) ff, w(1:nparm) f0(iv) = ff DO ii=1,nparm x0(ii,iv) = y(ii) 280 IF (n1 > 20) GO TO 470 DO ii=1,nparm x1(ii,n1) = x(ii,i1) xcl(ii,1) = x(ii,i1) (78 of 114) 오전 10:21:35

79 f1(n1) = f(i1) fcl(1) = f(i1) ic1(n1) = iv icj = iv GO TO 305 NEW LOCAL MINIMUM 290 nc = nc+1 ncp = ncp+1 WRITE(ipr,907) nc, ff, nfe1 WRITE(*,907) nc, ff, nfe1 907 FORMAT(' *** THE LOCAL MINIMUM NO. ', i2, ': ', g14.8, ', NFEV=', i5) DO ii=1,nparm w(ii) = y(ii)*mmax(ii) + mmin(ii) WRITE(ipr,903) ff, w(1:nparm) WRITE(*,903) ff, w(1:nparm) DO ii=1,nparm x0(ii,nc) = y(ii) xcl(ii,1) = y(ii) fcl(1) = ff f0(nc) = ff IF (nc >= 20) GO TO 475 it = 1 icj = nc CLUSTERING TO THE NEW POINT 305 nfe = nfe + nfe1 ic(i1) = icj ig = ig + 1 IF (ig >= ng) EXIT i = 1 in1 = i 310 DO j=1,ng IF (ic(j) /= 0) CYCLE IF (fcl(i) >= f(j)) CYCLE DO l1=1,nparm w(l1) = ABS(xcl(l1,i) - x(l1,j)) a = zero DO l1=1,nparm IF (w(l1) > a) a = w(l1) IF (a >= b) CYCLE in1 = in1 + 1 DO ii=1,nparm xcl(ii,in1) = x(ii,j) fcl(in1) = f(j) ic(j) = icj WRITE(ipr,902) icj WRITE(*,902) icj DO ii=1,nparm w(ii) = x(ii,j)*mmax(ii) + mmin(ii) WRITE(ipr,903) f(j), w(1:nparm) WRITE(*,903) f(j), w(1:nparm) ig = ig+1 IF (ig >= ng) EXIT (79 of 114) 오전 10:21:35

80 i = i+1 IF (i < in1) GO TO 310 IF (it /= 0) GO TO 20 PRINT RESULTS 395 WRITE(ipr,908) WRITE(*,908) 908 FORMAT(///' LOCAL MINIMA FOUND:'//) IF (nc <= 1) GO TO 430 inum = nc - 1 DO i=1,inum im = i fm = f0(i) inum1 = i + 1 DO j=inum1,nc IF (f0(j) >= fm) CYCLE im = j fm = f0(j) IF (im <= i) CYCLE a = fm y(1:nparm) = x0(1:nparm,im) inum1 = im-i DO j=1,inum1 inum2 = im - j f0(inum2+1) = f0(inum2) DO jj=1,nparm x0(jj,inum2+1) = x0(jj,inum2) f0(i) = a x0(1:nparm,i) = y(1:nparm) 430 IF (nc <= 0) GO TO 445 DO i=1,nc x0(1:nparm,i) = x0(1:nparm,i)*mmax(1:nparm) + mmin(1:nparm) WRITE(ipr,903) f0(i), x0(1:nparm,i) WRITE(*,903) f0(i), x0(1:nparm,i) 445 WRITE(ipr,911) nfe WRITE(*,911) nfe 911 FORMAT(//' NORMAL TERMINATION AFTER ', i5, ' FUNCTION EVALUATIONS'//) RETURN 460 WRITE(ipr,914) WRITE(*,914) 914 FORMAT(' *** DATA ERROR') STOP 465 WRITE(ipr,915) WRITE(*,915) 915 FORMAT(' *** TOO MANY SAMPLES') GO TO WRITE(ipr,916) WRITE(*,916) 916 FORMAT(' *** TOO MANY NEW SEED POINTS') GO TO (80 of 114) 오전 10:21:35

81 475 WRITE(ipr,917) WRITE(*,917) 917 FORMAT(' *** TOO MANY CLUSTERS') GO TO 395 END SUBROUTINE global SUBROUTINE fun(r, f, nparm, m, mmin, mmax) REAL (dp), INTENT(IN) :: r(:) REAL (dp), INTENT(OUT) :: f INTEGER, INTENT(IN) :: nparm INTEGER, INTENT(IN) :: m REAL (dp), INTENT(IN) :: mmin(:) REAL (dp), INTENT(IN) :: mmax(:) INTERFACE SUBROUTINE funct(x, f, nparm, m) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: x(:) REAL (dp), INTENT(OUT) :: f INTEGER, INTENT(IN) :: nparm, m END SUBROUTINE funct END INTERFACE REAL (dp) :: x(nparm) N.B. mmin = mid-point between lower & upper limits mmax = (upper - lower limits) / 2 x(1:nparm) = mmax(1:nparm)*r(1:nparm) + mmin(1:nparm) CALL funct(x, f, nparm, m) RETURN END SUBROUTINE fun ROUTINE NAME - LOCAL Code converted using TO_F90 by Alan Miller Date: Time: 13:16: LATEST REVISION - JULY 31, 1986 PURPOSE - MINIMUM OF A FUNCTION OF N VARIABLES USING A QUASI-NEWTON METHOD USAGE - CALL LOCAL (M, N, EPS, MAXFN, X, F, NFEV, mmin, mmax) ARGUMENTS M - THE NUMBER OF RESIDUAL FUNCTIONS (INPUT) NOT USED IN THIS ROUTINE. N - THE NUMBER OF PARAMETERS (I.E., THE LENGTH OF X) (INPUT) EPS - CONVERGENCE CRITERION. (INPUT). THE ACCURACY REQUIRED IN THE PARAMETER ESTIMATES. THIS CONVERGENCE CONDITION IS SATISFIED IF ON TWO SUCCESSIVE (81 of 114) 오전 10:21:35

82 ITERATIONS, THE PARAMETER ESTIMATES (I.E.,X(I), I=1,...,N) DIFFERS, COMPONENT BY COMPONENT, BY AT MOST EPS. MAXFN - MAXIMUM NUMBER OF FUNCTION EVALUATIONS (I.E., CALLS TO SUBROUTINE FUN) ALLOWED. (INPUT) X - VECTOR OF LENGTH N CONTAINING PARAMETER VALUES. ON INPUT, X MUST CONTAIN THE INITIAL PARAMETER ESTIMATES. ON OUTPUT, X CONTAINS THE FINAL PARAMETER ESTIMATES AS DETERMINED BY LOCAL. F - A SCALAR CONTAINING THE VALUE OF THE FUNCTION AT THE FINAL PARAMETER ESTIMATES. (OUTPUT) NFEV - THE NUMBER OF FUNCTION EVALUATIONS (OUTPUT) W - A VECTOR OF LENGTH 3*N USED AS WORKING SPACE. mmin - A VECTOR OF LENGTH N CONTAINING THE LOWER BOUNDS OF THE PARAMETERS, SO X(I) IS SEARCHED IN THE INTERVAL (mmin(i),mmax(i)). (INPUT) mmax - A VECTOR OF LENGTH N CONTAINING THE UPPER BOUNDS OF THE PARAMETERS. (INPUT) REQUIRED ROUTINES - UPDATE, FUN FUN - A USER SUPPLIED SUBROUTINE WHICH CALCULATES THE FUNCTION F FOR GIVEN PARAMETER VALUES X(1),X(2),...,X(N). THE CALLING SEQUENCE HAS THE FOLLOWING FORM CALL FUN(X, F, N, M, mmin, mmax) WHERE X IS A VEKTOR OF LENGTH N. FUN MUST APPEAR IN AN EXTERNAL STATEMENT IN THE CALLING PROGRAM. FUN MUST NOT ALTER THE VALUES OF X(I), I=1,...,N OR N SUBROUTINE local (m, n, eps, maxfn, x, f, nfev, mmin, mmax) N.B. Argument W has been removed. SPECIFICATIONS FOR ARGUMENTS INTEGER, INTENT(IN) :: m INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN) :: eps INTEGER, INTENT(IN) :: maxfn REAL (dp), INTENT(IN OUT) :: x(:) REAL (dp), INTENT(OUT) :: f INTEGER, INTENT(OUT) :: nfev REAL (dp), INTENT(IN) :: mmin(:) REAL (dp), INTENT(IN) :: mmax(:) SPECIFICATIONS FOR LOCAL VARIABLES REAL (dp) :: w(3*n) INTEGER :: ig, igg, im1, is, idiff, ir, ij, i, iopt, j, nm1, jj, jp1, l, & kj, k, link, itn, ii, jnt, np1, jb, nj, ier REAL (dp) :: hh, hjj, v, df, relx, gs0, diff, aeps, alpha, ff, tot, f1, f2, & z, gys, dgs, sig, zz, hhh, ghh, g(n), h(n*n) REAL (dp), PARAMETER :: reps = E-07_dp, ax = 0.1_dp, zero = 0.0_dp, & one = 1.0_dp, half = 0.5_dp, seven = 7.0_dp, & five = 5.0_dp, twelve = 12.0_dp, p1 = 0.1_dp INITIALIZATION (82 of 114) 오전 10:21:35

83 FIRST EXECUTABLE STATEMENT iopt = 0 IOPT - OPTIONS SELECTOR. (INPUT) IOPT = 0 CAUSES LOCAL TO INITIALIZE THE HESSIAN MATRIX H TO THE IDENTITY MATRIX. IOPT = 1 INDICATES THAT H HAS BEEN INITIALIZED BY THE USER TO A POSITIVE DEFINITE MATRIX. IOPT = 2 CAUSES LOCAL TO COMPUTE THE DIAGONAL VALUES OF THE HESSIAN MATRIX AND SET H TO A DIAGONAL MATRIX CONTAINING THESE VALUES. IOPT = 3 CAUSES LOCAL TO COMPUTE AN ESTIMATE OF THE HESSIAN IN H. ier = 0 hh = SQRT(reps) ig = n igg = n+n is = igg idiff = 1 ir = n w(1) = -one w(2) = zero w(3) = zero EVALUATE FUNCTION AT STARTING POINT g(1:n) = x(1:n) CALL fun (g, f, n, m, mmin, mmax) nfev = 1 IF (iopt == 1) GO TO 45 SET OFF-DIAGONAL ELEMENTS OF H TO 0.0 IF (n == 1) GO TO 25 ij = 2 DO i=2,n DO j=2,i h(ij) = zero ij = ij + 1 ij = ij + 1 IF (iopt /= 0) GO TO 25 SET DIAGONAL ELEMENTS OF H TO ONE ij = 0 DO i=1,n ij = ij+i h(ij) = one GO TO 80 GET DIAGONAL ELEMENTS OF HESSIAN 25 im1 = 1 nm1 = 1 np1 = n+1 DO i=2,np1 hhh = hh*max(abs(x(im1)), ax) g(im1) = x(im1) + hhh CALL fun (g, f2, n, m, mmin, mmax) g(im1) = g(im1) + hhh CALL fun (g, ff, n, m, mmin, mmax) h(nm1) = (ff-f2+f-f2)/(hhh*hhh) (83 of 114) 오전 10:21:35

84 g(im1) = x(im1) im1 = i nm1 = i+nm1 nfev = nfev+n+n IF (iopt /= 3.OR. n == 1) GO TO 45 GET THE REST OF THE HESSIAN jj = 1 ii = 2 DO i=2,n ghh = hh*max(abs(x(i)), ax) g(i) = x(i)+ghh CALL fun (g, f2, n, m, mmin, mmax) DO j=1,jj hhh = hh*max(abs(x(j)), ax) g(j) = x(j) + hhh CALL fun (g, ff, n, m, mmin, mmax) g(i) = x(i) CALL fun (g, f1, n, m, mmin, mmax) H(II) = (FF-F1-F2+F)*SQREPS h(ii) = (ff-f1-f2+f)/(hhh*ghh) ii = ii+1 g(j) = x(j) jj = jj+1 ii = ii+1 nfev = nfev + ((n*n-n)/2) FACTOR H TO L*D*L-TRANSPOSE 45 ir = n IF (n > 1) GO TO 50 IF (h(1) > zero) GO TO 80 h(1) = zero ir = 0 GO TO nm1 = n-1 jj = 0 DO j=1,n jp1 = j+1 jj = jj+j hjj = h(jj) IF (hjj > zero) GO TO 55 h(jj) = zero ir = ir-1 CYCLE 55 IF (j == n) CYCLE ij = jj l = 0 DO i=jp1,n l = l+1 ij = ij+i-1 v = h(ij)/hjj kj = ij DO k=i,n h(kj+l) = h(kj+l) - h(kj)*v kj = kj + k h(ij) = v (84 of 114) 오전 10:21:36

85 75 IF (ir /= n) THEN ier = 129 GO TO itn = 0 df = -one EVALUATE GRADIENT W(IG+I),I=1,...,N 85 link = 1 GO TO 260 BEGIN ITERATION LOOP 90 IF (nfev >= maxfn) GO TO 225 itn = itn+1 DO i=1,n w(i) = -w(ig+i) DETERMINE SEARCH DIRECTION W BY SOLVING H*W = -G WHERE H = L*D*L-TRANSPOSE IF (ir < n) GO TO 125 N.EQ. 1 g(1) = w(1) IF (n > 1) GO TO 100 w(1) = w(1)/h(1) GO TO 125 N > ii = 1 SOLVE L*W = -G DO i=2,n ij = ii ii = ii+i v = w(i) im1 = i-1 DO j=1,im1 ij = ij+1 v = v - h(ij)*w(j) g(i) = v w(i) = v SOLVE (D*LT)*Z = W WHERE LT = L-TRANSPOSE w(n) = w(n)/h(ii) jj = ii nm1 = n-1 DO nj=1,nm1 J = N-1,N-2,...,1 j = n-nj jp1 = j+1 jj = jj-jp1 v = w(j)/h(jj) ij = jj DO i=jp1,n ij = ij+i-1 v = v - h(ij)*w(i) w(j) = v (85 of 114) 오전 10:21:36

86 DETERMINE STEP LENGTH ALPHA 125 relx = zero gs0 = zero DO i=1,n w(is+i) = w(i) diff = ABS(w(i)) / MAX(ABS(x(i)),ax) relx = MAX(relx, diff) gs0 = gs0 + w(ig+i)*w(i) IF (relx == zero) GO TO 230 aeps = eps/relx ier = 130 IF (gs0 >= zero) GO TO 230 IF (df == zero) GO TO 230 ier = 0 alpha = (-df-df)/gs0 IF (alpha <= zero) alpha = one alpha = MIN(alpha,one) IF (idiff == 2) alpha = MAX(p1,alpha) ff = f tot = zero jnt = 0 SEARCH ALONG X + ALPHA*W 135 IF (nfev >= maxfn) GO TO 225 DO i=1,n w(i) = x(i) + alpha*w(is+i) CALL fun (w, f1, n, m, mmin, mmax) nfev = nfev+1 IF (f1 >= f) GO TO 165 f2 = f tot = tot + alpha 145 ier = 0 f = f1 x(1:n) = w(1:n) IF (jnt-1 < 0) THEN GO TO 155 ELSE IF (jnt-1 == 0) THEN GO TO 185 ELSE GO TO IF (nfev >= maxfn) GO TO 225 DO i=1,n w(i) = x(i) + alpha*w(is+i) CALL fun (w, f1, n, m, mmin, mmax) nfev = nfev+1 IF (f1 >= f) GO TO 190 IF (f1+f2 >= f+f.and. seven*f1+five*f2 > twelve*f) jnt = 2 tot = tot + alpha alpha = alpha + alpha GO TO IF (f == ff.and. idiff == 2.AND. relx > eps) ier = 130 IF (alpha < aeps) GO TO 230 IF (nfev >= maxfn) GO TO 225 alpha = half*alpha (86 of 114) 오전 10:21:36

87 DO i=1,n w(i) = x(i) + alpha*w(is+i) CALL fun (w, f2, n, m, mmin, mmax) nfev = nfev + 1 IF (f2 >= f) GO TO 180 tot = tot + alpha ier = 0 f = f2 x(1:n) = w(1:n) GO TO z = p1 IF (f1+f > f2+f2) z = one + half*(f-f1)/(f+f1-f2-f2) z = MAX(p1,z) alpha = z*alpha jnt = 1 GO TO IF (tot < aeps) GO TO alpha = tot SAVE OLD GRADIENT DO i=1,n w(i) = w(ig+i) EVALUATE GRADIENT W(IG+I), I=1,...,N link = 2 GO TO IF (nfev >= maxfn) GO TO 225 gys = zero DO i=1,n gys = gys + w(ig+i)*w(is+i) w(igg+i) = w(i) df = ff-f dgs = gys-gs0 IF (dgs <= zero) GO TO 90 IF (dgs + alpha*gs0 > zero) GO TO 215 UPDATE HESSIAN H USING COMPLEMENTARY DFP FORMULA sig = one/gs0 ir = -ir CALL update (h, n, w, sig, g, ir, 0, zero) DO i=1,n g(i) = w(ig+i) - w(igg+i) sig = one/(alpha*dgs) ir = -ir CALL update (h, n, g, sig, w, ir, 0, zero) GO TO 90 UPDATE HESSIAN USING DFP FORMULA 215 zz = alpha/(dgs - alpha*gs0) sig = -zz CALL update (h, n, w, sig, g, ir, 0, reps) z = dgs*zz - one DO i=1,n g(i) = w(ig+i) + z*w(igg+i) (87 of 114) 오전 10:21:36

88 sig = one/(zz*dgs*dgs) CALL update (h, n, g, sig, w, ir, 0, zero) GO TO 90 MAXFN FUNCTION EVALUATIONS 225 GO TO IF (idiff == 2) GO TO 235 CHANGE TO CENTRAL DIFFERENCES idiff = 2 GO TO IF (relx > eps.and. ier == 0) GO TO 85 COMPUTE H = L*D*L-TRANSPOSE AND OUTPUT IF (n == 1) GO TO 9000 np1 = n+1 nm1 = n-1 jj = (n*(np1))/2 DO jb=1,nm1 jp1 = np1-jb jj = jj-jp1 hjj = h(jj) ij = jj l = 0 DO i=jp1,n l = l+1 ij = ij+i-1 v = h(ij)*hjj kj = ij DO k=i,n h(kj+l) = h(kj+l) + h(kj)*v kj = kj + k h(ij) = v hjj = h(jj) GO TO 9000 EVALUATE GRADIENT 260 IF (idiff == 2) GO TO 270 FORWARD DIFFERENCES GRADIENT = W(IG+I), I=1,...,N DO i=1,n z = hh*max(abs(x(i)),ax) zz = x(i) x(i) = zz + z CALL fun (x, f1, n, m, mmin, mmax) w(ig+i) = (f1-f)/z x(i) = zz nfev = nfev+n SELECT CASE ( link ) CASE ( 1) GO TO 90 CASE ( 2) GO TO 200 END SELECT CENTRAL DIFFERENCES (88 of 114) 오전 10:21:36

89 GRADIENT = W(IG+I), I=1,...,N 270 DO i=1,n z = hh*max(abs(x(i)), ax) zz = x(i) x(i) = zz + z CALL fun (x, f1, n, m, mmin, mmax) x(i) = zz-z CALL fun (x, f2, n, m, mmin, mmax) w(ig+i) = (f1-f2)/(z+z) x(i) = zz nfev = nfev+n+n SELECT CASE ( link ) CASE ( 1) GO TO 90 CASE ( 2) GO TO 200 END SELECT RETURN 9000 RETURN END SUBROUTINE local ROUTINE NAME - UPDATE LATEST REVISION - JULY 31, 1986 (CHANGES IN COMMENTS) PURPOSE - NUCLEUS CALLED ONLY BY ROUTINE LOCAL REQD. ROUTINES - NONE REQUIRED SUBROUTINE update (a, n, z, sig, w, ir, mk, eps) SPECIFICATIONS FOR ARGUMENTS REAL (dp), INTENT(OUT) :: a(:) INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: z(:) REAL (dp), INTENT(IN) :: sig REAL (dp), INTENT(IN OUT) :: w(:) INTEGER, INTENT(OUT) :: ir INTEGER, INTENT(IN) :: mk REAL (dp), INTENT(IN) :: eps SPECIFICATIONS FOR LOCAL VARIABLES INTEGER :: j, jj, ij, jp1, i, ii, mm REAL (dp) :: ti, v, tim, al, r, b, gm, y REAL (dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp, four = 4.0_dp UPDATE FACTORS GIVEN IN A SIG*Z*Z-TRANSPOSE IS ADDED FIRST EXECUTABLE STATEMENT IF (n > 1) GO TO 5 (89 of 114) 오전 10:21:36

90 N.EQ. 1 a(1) = a(1) + sig*z(1)*z(1) ir = 1 IF (a(1) > zero) GO TO 9005 a(1) = zero ir = 0 GO TO 9005 N > 1 5 IF (sig > zero) GO TO 65 IF (sig == zero.or. ir == 0) GO TO 9005 ti = one/sig jj = 0 IF (mk == 0) GO TO 15 L*W = Z ON INPUT DO j=1,n jj = jj+j IF (a(jj) /= zero) ti = ti + (w(j)*w(j))/a(jj) GO TO 40 SOLVE L*W = Z 15 w(1:n) = z(1:n) DO j=1,n jj = jj+j v = w(j) IF (a(jj) > zero) GO TO 25 w(j) = zero CYCLE 25 ti = ti+(v*v)/a(jj) IF (j == n) CYCLE ij = jj jp1 = j+1 DO i=jp1,n ij = ij+i-1 w(i) = w(i) - v*a(ij) SET TI, TIM AND W 40 IF (ir <= 0) GO TO 45 IF (ti > zero) GO TO 50 IF (mk-1 > 0) THEN GO TO 55 ELSE GO TO ti = zero ir = -ir-1 GO TO ti = eps/sig IF (eps == zero) ir = ir-1 55 tim = ti ii = jj i = n DO j=1,n IF (a(ii) /= zero) tim = ti - (w(i)*w(i))/a(ii) w(i) = ti (90 of 114) 오전 10:21:36

91 ti = tim ii = ii-i i = i-1 mm = 1 GO TO mm = 0 tim = one/sig 70 jj = 0 UPDATE A DO j=1,n jj = jj+j ij = jj jp1 = j+1 UPDATE A(J,J) v = z(j) IF (a(jj) > zero) GO TO 85 A(J,J).EQ. ZERO IF (ir > 0.OR. sig < zero.or. v == zero) GO TO 80 ir = 1-ir a(jj) = (v*v)/tim IF (j == n) GO TO 9005 DO i=jp1,n ij = ij+i-1 a(ij) = z(i)/v GO TO ti = tim CYCLE A(J,J).GT. ZERO 85 al = v/a(jj) ti = w(j) IF (mm == 0) ti = tim + v*al r = ti/tim a(jj) = r*a(jj) IF (r == zero) EXIT IF (j == n) EXIT UPDATE REMAINDER OF COLUMN J b = al/ti IF (r > four) GO TO 95 DO i=jp1,n ij = ij+i-1 z(i) = z(i) - v*a(ij) a(ij) = a(ij) + b*z(i) GO TO gm = tim/ti DO i=jp1,n ij = ij+i-1 y = a(ij) a(ij) = b*z(i) + y*gm z(i) = z(i) - v*y 105 tim = ti IF (ir < 0) ir = -ir (91 of 114) 오전 10:21:36

92 9005 RETURN END SUBROUTINE update END MODULE global_minimum 예제 (5) vector From Mon Aug 14 14:28: Subject: vectors and 4-vectors Some vector and four-vector algebra. The possibility of defining vector entities with specific algebra makes the code more readable, compared to f77. A couple of examples: scalar product of two vectors: F90: s = v1*v2 F77: s = v1(1)*v2(1) + v1(2)*v2(2) + v1(3)*v2(3) vector product of two vectors: F90: V = A.X. B F77: v(1) = A(2)*B(3) - B(2)*A(3) V(2) = A(3)*B(1) - B(3)*A(1) V(3) = A(1)*B(2) - B(1)*A(2) this example is of interests for physics : mass of particle with four momentum p F90: m = sqrt(p*p) F77: s = p(1)**2 + p(2)**2 + p(3)**2 m = sqrt( p(4)**2 - s ) =========================================================================== =========================================================================== MODULE P3V_HANDLING type defined: P3V : vector with components (x, y, z) MODULE P4V_HANDLING type defined: P4V : 4-vector with components (E, V)=(E, V%x, V%y, V%z) metrix=diag(1,-1,-1,-1) First version: 19-AUG-1992 A. Bay [email protected] T. Monteiro Last update: 14-Aug-1995 =========================================================================== operations defined on vectors (v, v1, v2): v1 + v2 add two vectors v1 - v2 subtract (92 of 114) 오전 10:21:36

93 v1*v2 scalar product v1.x. v2 vector product v*r vector*real v = R 3 components of the vector set to same real value R v = R(3) vector equated to real(3) R(3) = v real(3) equated to vector functions: v2 = norm_3v(v1) vector v2 is v1 normalized to 1. ang = find_ang_3v(v1,v2) angle in [rad] between v1 and v2... operations on 4-vectors : v1 + v2 v1 - v2 v1*v2 scalar product v*r 4vector*real v = R 4 components of the 4vector set to same real value R v = R(4) vector equated to real(4) R(4) = v real(4) equated to vector functions: find_ang(v1,v2) angle in rad between v1%v and v2%v =========================================================================== MODULE P3V_HANDLING TYPE P3V REAL X REAL Y REAL Z END TYPE P3V INTERFACE ASSIGNMENT(=) MODULE PROCEDURE P3V_EQU_P3V, P3V_EQU_R, P3V_EQU_V3, V3_EQU_P3V END INTERFACE INTERFACE OPERATOR(+) MODULE PROCEDURE ADD_P3V END INTERFACE INTERFACE OPERATOR(-) MODULE PROCEDURE SUB_P3V END INTERFACE INTERFACE OPERATOR(*) MODULE PROCEDURE MULT_P3V,MULT_P3V_R,MULT_R_P3V END INTERFACE INTERFACE OPERATOR(.X.) MODULE PROCEDURE VEC_P3V END INTERFACE CONTAINS SUBROUTINE P3V_EQU_P3V(A,B) TYPE(P3V), INTENT(OUT) :: A TYPE(P3V), INTENT(IN) :: B A%X = B%X; A%Y = B%Y; A%Z = B%Z END SUBROUTINE P3V_EQU_P3V SUBROUTINE P3V_EQU_R(A,B) TYPE(P3V), INTENT(OUT) :: A REAL, INTENT(IN) :: B A%X = B; A%Y = B; A%Z = B END SUBROUTINE P3V_EQU_R SUBROUTINE P3V_EQU_V3(A,R) (93 of 114) 오전 10:21:36

94 TYPE(P3V), INTENT(OUT):: A REAL, DIMENSION(3), INTENT(IN) :: R A%X = R(1); A%Y = R(2); A%Z = R(3) END SUBROUTINE P3V_EQU_V3 SUBROUTINE V3_EQU_P3V(R,A) TYPE(P3V), INTENT(IN):: A REAL, DIMENSION(3), INTENT(OUT) :: R R(1) = A%X; R(2) = A%Y; R(3) = A%Z END SUBROUTINE V3_EQU_P3V FUNCTION ADD_P3V(A,B) TYPE(P3V) ADD_P3V TYPE(P3V), INTENT(IN):: A,B ADD_P3V%X = A%X + B%X; ADD_P3V%Y = A%Y + B%Y; ADD_P3V%Z = A%Z + B%Z END FUNCTION ADD_P3V FUNCTION SUB_P3V(A,B) TYPE(P3V) SUB_P3V TYPE(P3V), INTENT(IN):: A,B SUB_P3V%X = A%X - B%X; SUB_P3V%Y = A%Y - B%Y; SUB_P3V%Z = A%Z - B%Z END FUNCTION SUB_P3V FUNCTION MULT_P3V(A,B) REAL MULT_P3V TYPE(P3V), INTENT(IN):: A,B MULT_P3V = A%X*B%X + A%Y*B%Y + A%Z*B%Z END FUNCTION MULT_P3V FUNCTION MULT_P3V_R(A,R) TYPE(P3V) MULT_P3V_R TYPE(P3V), INTENT(IN) :: A REAL, INTENT(IN) :: R MULT_P3V_R%X = A%X * R; MULT_P3V_R%Y = A%Y * R; MULT_P3V_R%Z = A%Z * R END FUNCTION MULT_P3V_R FUNCTION MULT_R_P3V(R,A) TYPE(P3V) MULT_R_P3V TYPE(P3V), INTENT(IN) :: A REAL, INTENT(IN) :: R MULT_R_P3V%X = A%X * R; MULT_R_P3V%Y = A%Y * R; MULT_R_P3V%Z = A%Z * R END FUNCTION MULT_R_P3V FUNCTION VEC_P3V(A,B) TYPE(P3V) VEC_P3V TYPE(P3V), INTENT(IN) :: A,B VEC_P3V%X=A%Y*B%Z-B%Y*A%Z VEC_P3V%Y=A%Z*B%X-B%Z*A%X VEC_P3V%Z=A%X*B%Y-B%X*A%Y END FUNCTION VEC_P3V FUNCTION NORM_3V(V) TYPE (P3V) ::NORM_3V TYPE (P3V),INTENT(IN) ::V REAL :: N N = SQRT(V*V) IF (N>0) THEN NORM_3V%X = V%X/N NORM_3V%Y = V%Y/N NORM_3V%Z = V%Z/N ENDIF END FUNCTION NORM_3V Angles FUNCTION FIND_ANG_3V(A,B) REAL ::FIND_ANG_3V TYPE (P3V), INTENT(IN) ::A,B (94 of 114) 오전 10:21:36

95 REAL ::AB,D,X AB = A*B; D = (A*A)*(B*B) IF (D==0.) THEN FIND_ANG_3V = ELSE X = (AB/SQRT(D)) IF (ABS(X)>1.) THEN WRITE(*,*)'Error in Find_ang', X WRITE(*,*)'A=',A WRITE(*,*)'B=',B IF (X<-1.) X = -1. IF (X>1.) X = 1. ENDIF FIND_ANG_3V = ACOS(X) ENDIF END FUNCTION FIND_ANG_3V END MODULE P3V_HANDLING ============================================================= MODULE P4V_HANDLING USE P3V_HANDLING TYPE P4V TYPE(P3V) V REAL E END TYPE P4V INTERFACE ASSIGNMENT(=) MODULE PROCEDURE P4V_EQU_P4V, P4V_EQU_R,P4V_EQU_V4,V4_EQU_P4V END INTERFACE INTERFACE OPERATOR(+) MODULE PROCEDURE ADD_P4V END INTERFACE INTERFACE OPERATOR(-) MODULE PROCEDURE SUB_P4V END INTERFACE INTERFACE OPERATOR(*) MODULE PROCEDURE MULT_P4V,MULT_P4V_R,MULT_R_P4V END INTERFACE INTERFACE FIND_ANG MODULE PROCEDURE FIND_ANG_3V,FIND_ANG_4V END INTERFACE CONTAINS SUBROUTINE P4V_EQU_P4V(A,B) TYPE(P4V), INTENT(IN) :: B TYPE(P4V), INTENT(OUT) :: A A%V = B%V; A%E = B%E END SUBROUTINE P4V_EQU_P4V SUBROUTINE P4V_EQU_R(A,B) TYPE(P4V), INTENT(OUT) :: A REAL, INTENT(IN) :: B A%V%X = B; A%V%Y = B; A%V%Z = B A%E = B END SUBROUTINE P4V_EQU_R SUBROUTINE P4V_EQU_V4(A,R) TYPE(P4V), INTENT(OUT):: A REAL, DIMENSION(4), INTENT(IN) :: R(4) A%V%X = R(1); A%V%Y = R(2); A%V%Z = R(3) (95 of 114) 오전 10:21:36

96 A%E = R(4) END SUBROUTINE P4V_EQU_V4 SUBROUTINE V4_EQU_P4V(R,A) TYPE(P4V), INTENT(IN) :: A REAL, DIMENSION(4), INTENT(OUT) :: R R(1) = A%V%X; R(2) = A%V%Y; R(3) = A%V%Z R(4) = A%E END SUBROUTINE V4_EQU_P4V FUNCTION ADD_P4V(A,B) TYPE(P4V) ADD_P4V TYPE(P4V), INTENT(IN) :: A,B ADD_P4V%V = A%V + B%V ADD_P4V%E = A%E + B%E END FUNCTION ADD_P4V FUNCTION SUB_P4V(A,B) TYPE(P4V) SUB_P4V TYPE(P4V), INTENT(IN) :: A,B SUB_P4V%V = A%V - B%V SUB_P4V%E = A%E - B%E END FUNCTION SUB_P4V FUNCTION MULT_P4V(A,B) REAL MULT_P4V TYPE(P4V), INTENT(IN) :: A,B MULT_P4V = - (A%V*B%V) + A%E*B%E END FUNCTION MULT_P4V FUNCTION MULT_P4V_R(A,R) TYPE(P4V) MULT_P4V_R TYPE(P4V), INTENT(IN) :: A REAL, INTENT(IN) :: R MULT_P4V_R%V = A%V * R; MULT_P4V_R%E = A%E * R END FUNCTION MULT_P4V_R FUNCTION MULT_R_P4V(R,A) TYPE(P4V) MULT_R_P4V TYPE(P4V), INTENT(IN) :: A REAL, INTENT(IN) :: R MULT_R_P4V%V = A%V * R; MULT_R_P4V%E = A%E * R END FUNCTION MULT_R_P4V angles FUNCTION FIND_ANG_4V(A,B) REAL ::FIND_ANG_4V TYPE (P4V), INTENT(IN) ::A,B FIND_ANG_4V = FIND_ANG_3V(A%V,B%V) END FUNCTION FIND_ANG_4V END MODULE P4V_HANDLING =================================================================================== =================================================================================== 예제 (1) simulated annealing program (96 of 114) 오전 10:21:36

97 This file is an example of the Corana et al. simulated annealing algorithm for multimodal and robust optimization as implemented and modified by Goffe, Ferrier and Rogers. Counting the above line ABSTRACT as 1, the routine itself (SA), with its supplementary routines, is on lines A multimodal example from Judge et al. (FCN) is on lines The rest of this file (lines 1-149) is a driver routine with values appropriate for the Judge example. Thus, this example is ready to run. To understand the algorithm, the documentation for SA on lines should be read along with the parts of the paper that describe simulated annealing. Then the following lines will then aid the user in becomming proficient with this implementation of simulated annealing. Learning to use SA: Use the sample function from Judge with the following suggestions to get a feel for how SA works. When you've done this, you should be ready to use it on most any function with a fair amount of expertise. 1. Run the program as is to make sure it runs okay. Take a look at the intermediate output and see how it optimizes as temperature (T) falls. Notice how the optimal point is reached and how falling T reduces VM. 2. Look through the documentation to SA so the following makes a bit of sense. In line with the paper, it shouldn't be that hard to figure out. The core of the algorithm is described on pp and on pp Also see Corana et al. pp To see how it selects points and makes decisions about uphill and downhill moves, set IPRINT = 3 (very detailed intermediate output) and MAXEVL = 100 (only 100 function evaluations to limit output). 4. To see the importance of different temperatures, try starting with a very low one (say T = 10E-5). You'll see (i) it never escapes from the local optima (in annealing terminology, it quenches) & (ii) the step length (VM) will be quite small. This is a key part of the algorithm: as temperature (T) falls, step length does too. In a minor point here, note how VM is quickly reset from its initial value. Thus, the input VM is not very important. This is all the more reason to examine VM once the algorithm is underway. 5. To see the effect of different parameters and their effect on the speed of the algorithm, try RT =.95 & RT =.1. Notice the vastly different speed for optimization. Also try NT = 20. Note that this sample function is quite easy to optimize, so it will tolerate big changes in these parameters. RT and NT are the parameters one should adjust to modify the runtime of the algorithm and its robustness. 6. Try constraining the algorithm with either LB or UB. MODULE simulated_anneal ABSTRACT: Simulated annealing is a global optimization method that distinguishes (97 of 114) 오전 10:21:36

98 between different local optima. Starting from an initial point, the algorithm takes a step and the function is evaluated. When minimizing a function, any downhill step is accepted and the process repeats from this new point. An uphill step may be accepted. Thus, it can escape from local optima. This uphill decision is made by the Metropolis criteria. As the optimization process proceeds, the length of the steps decline and the algorithm closes in on the global optimum. Since the algorithm makes very few assumptions regarding the function to be optimized, it is quite robust with respect to non-quadratic surfaces. The degree of robustness can be adjusted by the user. In fact, simulated annealing can be used as a local optimizer for difficult functions. This implementation of simulated annealing was used in "Global Optimizatio of Statistical Functions with Simulated Annealing," Goffe, Ferrier and Rogers, Journal of Econometrics, vol. 60, no. 1/2, Jan./Feb. 1994, pp Briefly, we found it competitive, if not superior, to multiple restarts of conventional optimization routines for difficult optimization problems. For more information on this routine, contact its author: Bill Goffe, [email protected] This version in Fortran 90 has been prepared by Alan Miller. It is compatible with Lahey's ELF90 compiler. N.B. The 3 last arguments have been removed from subroutine sa. these were work arrays and are now internal to the routine. bigpond.net.au URL : Latest revision of Fortran 90 version - 3 August 1997 IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) The following variables were in COMMON /raset1/ REAL, SAVE :: u(97), cc, cd, cm INTEGER, SAVE :: i97, j97 CONTAINS SUBROUTINE sa(n, x, max, rt, eps, ns, nt, neps, maxevl, lb, ub, c, iprint, & iseed1, iseed2, t, vm, xopt, fopt, nacc, nfcnev, nobds, ier) Version: 3.2 Date: 1/22/94. (98 of 114) 오전 10:21:36

99 Differences compared to Version 2.0: 1. If a trial is out of bounds, a point is randomly selected from LB(i) to UB(i). Unlike in version 2.0, this trial is evaluated and is counted in acceptances and rejections. All corresponding documentation was changed as well. Differences compared to Version 3.0: 1. If VM(i) > (UB(i) - LB(i)), VM is set to UB(i) - LB(i). The idea is that if T is high relative to LB & UB, most points will be accepted, causing VM to rise. But, in this situation, VM has little meaning; particularly if VM is larger than the acceptable region. Setting VM to this size still allows all parts of the allowable region to be selected. Differences compared to Version 3.1: 1. Test made to see if the initial temperature is positive. 2. WRITE statements prettied up. 3. References to paper updated. Synopsis: This routine implements the continuous simulated annealing global optimization algorithm described in Corana et al.'s article "Minimizing Multimodal Functions of Continuous Variables with the "Simulated Annealing" Algorithm" in the September 1987 (vol. 13, no. 3, pp ) issue of the ACM Transactions on Mathematical Software. A very quick (perhaps too quick) overview of SA: SA tries to find the global optimum of an N dimensional function. It moves both up and downhill and as the optimization process proceeds, it focuses on the most promising area. To start, it randomly chooses a trial point within the step length VM (a vector of length N) of the user selected starting point. The function is evaluated at this trial point and its value is compared to its value at the initial point. In a maximization problem, all uphill moves are accepted and the algorithm continues from that trial point. Downhill moves may be accepted; the decision is made by the Metropolis criteria. It uses T (temperature) and the size of the downhill move in a probabilistic manner. The smaller T and the size of the downhill move are, the more likely that move will be accepted. If the trial is accepted, the algorithm moves on from that point. If it is rejected, another point is chosen instead for a trial evaluation. Each element of VM periodically adjusted so that half of all function evaluations in that direction are accepted. A fall in T is imposed upon the system with the RT variable by T(i+1) = RT*T(i) where i is the ith iteration. Thus, as T declines, downhill moves are less likely to be accepted and the percentage of rejections rise. Given the scheme for the selection for VM, VM falls. Thus, as T declines, VM falls and SA focuses upon the most promising area for optimization. The importance of the parameter T: The parameter T is crucial in using SA successfully. It influences VM, the step length over which the algorithm searches for optima. For a small intial T, the step length may be too small; thus not enough of the function might be evaluated to find the global optima. The user should carefully examine VM in the intermediate output (set IPRINT = 1) to make sure that VM is appropriate. The relationship between the initial temperature and the resulting step length is function dependent. (99 of 114) 오전 10:21:36

100 To determine the starting temperature that is consistent with optimizing a function, it is worthwhile to run a trial run first. Set RT = 1.5 and T = 1.0. With RT > 1.0, the temperature increases and VM rises as well. Then select the T that produces a large enough VM. For modifications to the algorithm and many details on its use, (particularly for econometric applications) see Goffe, Ferrier and Rogers, "Global Optimization of Statistical Functions with Simulated Annealing," Journal of Econometrics, vol. 60, no. 1/2, Jan./Feb. 1994, pp For more information, contact Bill Goffe Department of Economics and International Business University of Southern Mississippi Hattiesburg, MS (601) (office) (601) (fax) [email protected] (Internet) As far as possible, the parameters here have the same name as in the description of the algorithm on pp of Corana et al. In this description, SP is single precision, DP is double precision, INT is integer, L is logical and (N) denotes an array of length n. Thus, DP(N) denotes a double precision array of length n. Input Parameters: Note: The suggested values generally come from Corana et al. To drastically reduce runtime, see Goffe et al., pp for suggestions on choosing the appropriate RT and NT. N - Number of variables in the function to be optimized. (INT) X - The starting values for the variables of the function to be optimized. (DP(N)) MAX - Denotes whether the function should be maximized or minimized. A true value denotes maximization while a false value denotes minimization. Intermediate output (see IPRINT) takes this into account. (L) RT - The temperature reduction factor. The value suggested by Corana et al. is.85. See Goffe et al. for more advice. (DP) EPS - Error tolerance for termination. If the final function values from the last neps temperatures differ from the corresponding value at the current temperature by less than EPS and the final function value at the current temperature differs from the current optimal function value by less than EPS, execution terminates and IER = 0 is returned. (EP) NS - Number of cycles. After NS*N function evaluations, each element of VM is adjusted so that approximately half of all function evaluations are accepted. The suggested value is 20. (INT) NT - Number of iterations before temperature reduction. After NT*NS*N function evaluations, temperature (T) is changed by the factor RT. Value suggested by Corana et al. is MAX(100, 5*N). See Goffe et al. for further advice. (INT) NEPS - Number of final function values used to decide upon termi- nation. See EPS. Suggested value is 4. (INT) MAXEVL - The maximum number of function evaluations. If it is exceeded, IER = 1. (INT) (100 of 114) 오전 10:21:36

101 LB - The lower bound for the allowable solution variables. (DP(N)) UB - The upper bound for the allowable solution variables. (DP(N)) If the algorithm chooses X(I).LT. LB(I) or X(I).GT. UB(I), I = 1, N, a point is from inside is randomly selected. This This focuses the algorithm on the region inside UB and LB. Unless the user wishes to concentrate the search to a particular region, UB and LB should be set to very large positive and negative values, respectively. Note that the starting vector X should be inside this region. Also note that LB and UB are fixed in position, while VM is centered on the last accepted trial set of variables that optimizes the function. C - Vector that controls the step length adjustment. The suggested value for all elements is 2.0. (DP(N)) IPRINT - controls printing inside SA. (INT) Values: 0 - Nothing printed. 1 - Function value for the starting value and summary results before each temperature reduction. This includes the optimal function value found so far, the total number of moves (broken up into uphill, downhill, accepted and rejected), the number of out of bounds trials, the number of new optima found at this temperature, the current optimal X and the step length VM. Note that there are N*NS*NT function evalutations before each temperature reduction. Finally, notice is is also given upon achieveing the termination criteria. 2 - Each new step length (VM), the current optimal X (XOPT) and the current trial X (X). This gives the user some idea about how far X strays from XOPT as well as how VM is adapting to the function. 3 - Each function evaluation, its acceptance or rejection and new optima. For many problems, this option will likely require a small tree if hard copy is used. This option is best used to learn about the algorithm. A small value for MAXEVL is thus recommended when using IPRINT = 3. Suggested value: 1 Note: For a given value of IPRINT, the lower valued options (other than 0) are utilized. ISEED1 - The first seed for the random number generator RANMAR. 0 <= ISEED1 <= (INT) ISEED2 - The second seed for the random number generator RANMAR. 0 <= ISEED2 <= Different values for ISEED1 and ISEED2 will lead to an entirely different sequence of trial points and decisions on downhill moves (when maximizing). See Goffe et al. on how this can be used to test the results of SA. (INT) Input/Output Parameters: T - On input, the initial temperature. See Goffe et al. for advice. On output, the final temperature. (DP) VM - The step length vector. On input it should encompass the region of interest given the starting value X. For point X(I), the next trial point is selected is from X(I) - VM(I) to X(I) + VM(I). (101 of 114) 오전 10:21:36

102 Since VM is adjusted so that about half of all points are accepted, the input value is not very important (i.e. is the value is off, SA adjusts VM to the correct value). (DP(N)) Output Parameters: XOPT - The variables that optimize the function. (DP(N)) FOPT - The optimal value of the function. (DP) NACC - The number of accepted function evaluations. (INT) NFCNEV - The total number of function evaluations. In a minor point, note that the first evaluation is not used in the core of the algorithm; it simply initializes the algorithm. (INT). NOBDS - The total number of trial function evaluations that would have been out of bounds of LB and UB. Note that a trial point is randomly selected between LB and UB. (INT) IER - The error return number. (INT) Values: 0 - Normal return; termination criteria achieved. 1 - Number of function evaluations (NFCNEV) is greater than the maximum number (MAXEVL). 2 - The starting value (X) is not inside the bounds (LB and UB). 3 - The initial temperature is not positive Should not be seen; only used internally. Work arrays that must be dimensioned in the calling routine: RWK1 (DP(NEPS)) (FSTAR in SA) RWK2 (DP(N)) (XP " " ) IWK (INT(N)) (NACP " " ) N.B. In the Fortran 90 version, these are automatic arrays. Required Functions (included): EXPREP - Replaces the function EXP to avoid under- and overflows. It may have to be modified for non IBM-type main- frames. (DP) RMARIN - Initializes the random number generator RANMAR. RANMAR - The actual random number generator. Note that RMARIN must run first (SA does this). It produces uniform random numbers on [0,1]. These routines are from Usenet's comp.lang.fortran. For a reference, see "Toward a Universal Random Number Generator" by George Marsaglia and Arif Zaman, Florida State University Report: FSU-SCRI (1987). It was later modified by F. James and published in "A Review of Pseudo-random Number Generators." For further information, contact [email protected]. These routines are designed to be portable on any machine with a 24-bit or more mantissa. I have found it produces identical results on a IBM 3081 and a Cray Y-MP. Required Subroutines (included): PRTVEC - Prints vectors. PRT1... PRT10 - Prints intermediate output. FCN - Function to be optimized. The form is SUBROUTINE FCN(N, X, F) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) INTEGER, INTENT(IN) :: N REAL (dp), INTENT(IN) :: X(N) REAL (dp), INTENT(OUT) :: F (102 of 114) 오전 10:21:36

103 ... function code with F = F(X)... RETURN END Note: This is the same form used in the multivariable minimization algorithms in the IMSL edition 10 library. Machine Specific Features: 1. EXPREP may have to be modified if used on non-ibm type main- frames. Watch for under- and overflows in EXPREP. 2. Some FORMAT statements use G25.18; this may be excessive for some machines. 3. RMARIN and RANMAR are designed to be protable; they should not cause any problems. REAL (dp), INTENT(IN) :: lb(:), ub(:), c(:), eps, rt REAL (dp), INTENT(IN OUT) :: x(:), t, vm(:) REAL (dp), INTENT(OUT) :: xopt(:), fopt INTEGER, INTENT(IN) :: n, ns, nt, neps, maxevl, iprint, iseed1, iseed2 INTEGER, INTENT(OUT) :: nacc, nfcnev, nobds, ier LOGICAL, INTENT(IN) :: max Type all internal variables. REAL (dp) :: f, fp, p, pp, ratio, xp(n), fstar(neps) INTEGER :: nup, ndown, nrej, nnew, lnobds, h, i, j, m, nacp(n) LOGICAL :: quit INTERFACE SUBROUTINE fcn(n, theta, h) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN) :: theta(:) REAL (dp), INTENT(OUT) :: h END SUBROUTINE fcn END INTERFACE Initialize the random number generator RANMAR. CALL rmarin(iseed1, iseed2) Set initial values. nacc = 0 nobds = 0 nfcnev = 0 ier = 99 DO i = 1, n xopt(i) = x(i) nacp(i) = 0 fstar = 1.0D+20 If the initial temperature is not positive, notify the user and return to the calling routine. IF (t <= 0.0) THEN WRITE(*,'(/, " THE INITIAL TEMPERATURE IS NOT POSITIVE. "/, & & " reset the variable t. "/)') ier = 3 RETURN (103 of 114) 오전 10:21:36

104 If the initial value is out of bounds, notify the user and return to the calling routine. DO i = 1, n IF ((x(i) > ub(i)).or. (x(i) < lb(i))) THEN CALL prt1() ier = 2 RETURN Evaluate the function with input X and return value as F. CALL fcn(n, x, f) If the function is to be minimized, switch the sign of the function. Note that all intermediate and final output switches the sign back to eliminate any possible confusion for the user. IF(.NOT. max) f = -f nfcnev = nfcnev + 1 fopt = f fstar(1) = f IF(iprint >= 1) CALL prt2(max, n, x, f) Start the main loop. Note that it terminates if (i) the algorithm succesfully optimizes the function or (ii) there are too many function evaluations (more than MAXEVL). 100 nup = 0 nrej = 0 nnew = 0 ndown = 0 lnobds = 0 DO m = 1, nt DO j = 1, ns DO h = 1, n Generate XP, the trial value of X. Note use of VM to choose XP. DO i = 1, n IF (i == h) THEN xp(i) = x(i) + (ranmar()* ) * vm(i) ELSE xp(i) = x(i) If XP is out of bounds, select a point in bounds for the trial. IF((xp(i) < lb(i)).or. (xp(i) > ub(i))) THEN xp(i) = lb(i) + (ub(i) - lb(i))*ranmar() lnobds = lnobds + 1 nobds = nobds + 1 IF(iprint >= 3) CALL prt3(max, n, xp, x, f) Evaluate the function with the trial point XP and return as FP. CALL fcn(n, xp, fp) IF(.NOT. max) fp = -fp nfcnev = nfcnev + 1 IF(iprint >= 3) CALL prt4(max, n, xp, x, fp, f) If too many function evaluations occur, terminate the algorithm. IF(nfcnev >= maxevl) THEN CALL prt5() (104 of 114) 오전 10:21:36

105 IF (.NOT. max) fopt = -fopt ier = 1 RETURN Accept the new point if the function value increases. IF(fp >= f) THEN IF(iprint >= 3) THEN WRITE(*,'(" POINT ACCEPTED")') x(1:n) = xp(1:n) f = fp nacc = nacc + 1 nacp(h) = nacp(h) + 1 nup = nup + 1 If greater than any other point, record as new optimum. IF (fp > fopt) THEN IF(iprint >= 3) THEN WRITE(*,'(" NEW OPTIMUM")') xopt(1:n) = xp(1:n) fopt = fp nnew = nnew + 1 If the point is lower, use the Metropolis criteria to decide on acceptance or rejection. ELSE p = exprep((fp - f)/t) pp = ranmar() IF (pp < p) THEN IF(iprint >= 3) CALL prt6(max) x(1:n) = xp(1:n) f = fp nacc = nacc + 1 nacp(h) = nacp(h) + 1 ndown = ndown + 1 ELSE nrej = nrej + 1 IF(iprint >= 3) CALL prt7(max) Adjust VM so that approximately half of all evaluations are accepted. DO i = 1, n ratio = DBLE(nacp(i)) /DBLE(ns) IF (ratio >.6) THEN vm(i) = vm(i)*(1. + c(i)*(ratio -.6)/.4) ELSE IF (ratio <.4) THEN vm(i) = vm(i)/(1. + c(i)*((.4 - ratio)/.4)) IF (vm(i) > (ub(i)-lb(i))) THEN vm(i) = ub(i) - lb(i) IF(iprint >= 2) THEN (105 of 114) 오전 10:21:36

106 CALL prt8(n, vm, xopt, x) nacp(1:n) = 0 IF(iprint >= 1) THEN CALL prt9(max,n,t,xopt,vm,fopt,nup,ndown,nrej,lnobds,nnew) Check termination criteria. quit =.false. fstar(1) = f IF ((fopt - fstar(1)) <= eps) quit =.true. DO i = 1, neps IF (ABS(f - fstar(i)) > eps) quit =.false. Terminate SA if appropriate. IF (quit) THEN x(1:n) = xopt(1:n) ier = 0 IF (.NOT. max) fopt = -fopt IF(iprint >= 1) CALL prt10() RETURN If termination criteria is not met, prepare for another loop. t = rt*t DO i = neps, 2, -1 fstar(i) = fstar(i-1) f = fopt x(1:n) = xopt(1:n) Loop again. GO TO 100 END SUBROUTINE sa FUNCTION exprep(rdum) RESULT(fn_val) This function replaces exp to avoid under- and overflows and is designed for IBM 370 type machines. It may be necessary to modify it for other machines. Note that the maximum and minimum values of EXPREP are such that they has no effect on the algorithm. REAL (dp), INTENT(IN) :: rdum REAL (dp) :: fn_val IF (rdum > 174._dp) THEN fn_val = 3.69D+75 ELSE IF (rdum < dp) THEN fn_val = 0.0_dp ELSE fn_val = EXP(rdum) RETURN (106 of 114) 오전 10:21:36

107 END FUNCTION exprep SUBROUTINE rmarin(ij, kl) This subroutine and the next function generate random numbers. See the comments for SA for more information. The only changes from the orginal code is that (1) the test to make sure that RMARIN runs first was taken out since SA assures that this is done (this test didn't compile under IBM's VS Fortran) and (2) typing ivec as integer was taken out since ivec isn't used. With these exceptions, all following lines are original. This is the initialization routine for the random number generator RANMAR() NOTE: The seed variables can have values between: 0 <= IJ <= <= KL <= INTEGER, INTENT(IN) :: ij, kl INTEGER :: i, j, k, l, ii, jj, m REAL :: s, t IF( ij < 0.OR. ij > OR. kl < 0.OR. kl > ) THEN WRITE(*, '(A)') ' The first random number seed must have a value ', & 'between 0 AND 31328' WRITE(*, '(A)') ' The second seed must have a value between 0 and 30081' STOP i = MOD(ij/177, 177) + 2 j = MOD(ij, 177) + 2 k = MOD(kl/169, 178) + 1 l = MOD(kl, 169) DO ii = 1, 97 s = 0.0 t = 0.5 DO jj = 1, 24 m = MOD(MOD(i*j, 179)*k, 179) i = j j = k k = m l = MOD(53*l+1, 169) IF (MOD(l*m, 64) >= 32) THEN s = s + t t = 0.5 * t u(ii) = s cc = / cd = / cm = / i97 = 97 j97 = 33 RETURN END SUBROUTINE rmarin FUNCTION ranmar() RESULT(fn_val) REAL :: fn_val Local variable REAL :: uni (107 of 114) 오전 10:21:36

108 uni = u(i97) - u(j97) IF( uni < 0.0 ) uni = uni u(i97) = uni i97 = i97-1 IF(i97 == 0) i97 = 97 j97 = j97-1 IF(j97 == 0) j97 = 97 cc = cc - cd IF( cc < 0.0 ) cc = cc + cm uni = uni - cc IF( uni < 0.0 ) uni = uni fn_val = uni RETURN END FUNCTION ranmar SUBROUTINE prt1() This subroutine prints intermediate output, as does PRT2 through PRT10. Note that if SA is minimizing the function, the sign of the function value and the directions (up/down) are reversed in all output to correspond with the actual function optimization. This correction is because SA was written to maximize functions and it minimizes by maximizing the negative a function. WRITE(*, '(/, " THE STARTING VALUE (X) IS OUTSIDE THE BOUNDS "/, & & " (lb AND ub). execution terminated without any"/, & & " optimization. respecify x, ub OR lb so that "/, & & " lb(i) < x(i) < ub(i), i = 1, n. "/)') RETURN END SUBROUTINE prt1 SUBROUTINE prt2(max, n, x, f) REAL (dp), INTENT(IN) :: x(:), f INTEGER, INTENT(IN) :: n LOGICAL, INTENT(IN) :: max WRITE(*, '(" ")') CALL prtvec(x,n,'initial X') IF (max) THEN WRITE(*, '(" INITIAL F: ",/, G25.18)') f ELSE WRITE(*, '(" INITIAL F: ",/, G25.18)') -f RETURN END SUBROUTINE prt2 SUBROUTINE prt3(max, n, xp, x, f) REAL (dp), INTENT(IN) :: xp(:), x(:), f INTEGER, INTENT(IN) :: n LOGICAL, INTENT(IN) :: max WRITE(*, '(" ")') CALL prtvec(x, n, 'CURRENT X') (108 of 114) 오전 10:21:36

109 IF (max) THEN WRITE(*, '(" CURRENT F: ", G25.18)') f ELSE WRITE(*, '(" CURRENT F: ", G25.18)') -f CALL prtvec(xp, n, 'TRIAL X') WRITE(*, '(" POINT REJECTED SINCE OUT OF BOUNDS")') RETURN END SUBROUTINE prt3 SUBROUTINE prt4(max, n, xp, x, fp, f) REAL (dp), INTENT(IN) :: xp(:), x(:), fp, f INTEGER, INTENT(IN) :: n LOGICAL, INTENT(IN) :: max WRITE(*,'(" ")') CALL prtvec(x,n,'current X') IF (max) THEN WRITE(*,'(" CURRENT F: ",G25.18)') f CALL prtvec(xp,n,'trial X') WRITE(*,'(" RESULTING F: ",G25.18)') fp ELSE WRITE(*,'(" CURRENT F: ",G25.18)') -f CALL prtvec(xp,n,'trial X') WRITE(*,'(" RESULTING F: ",G25.18)') -fp RETURN END SUBROUTINE prt4 SUBROUTINE prt5() WRITE(*, '(/, " TOO MANY FUNCTION EVALUATIONS; CONSIDER "/, & & " increasing maxevl OR eps, OR decreasing "/, & & " nt OR rt. these results are likely TO be "/, " poor.",/)') RETURN END SUBROUTINE prt5 SUBROUTINE prt6(max) LOGICAL, INTENT(IN) :: max IF (max) THEN WRITE(*,'(" THOUGH LOWER, POINT ACCEPTED")') ELSE WRITE(*,'(" THOUGH HIGHER, POINT ACCEPTED")') RETURN END SUBROUTINE prt6 SUBROUTINE prt7(max) LOGICAL, INTENT(IN) :: max (109 of 114) 오전 10:21:36

110 IF (max) THEN WRITE(*,'(" LOWER POINT REJECTED")') ELSE WRITE(*,'(" HIGHER POINT REJECTED")') RETURN END SUBROUTINE prt7 SUBROUTINE prt8(n, vm, xopt, x) REAL (dp), INTENT(IN) :: vm(:), xopt(:), x(:) INTEGER, INTENT(IN) :: n WRITE(*,'(/, " intermediate results after step length adjustment", /)') CALL prtvec(vm, n, 'NEW STEP LENGTH (VM)') CALL prtvec(xopt, n, 'CURRENT OPTIMAL X') CALL prtvec(x, n, 'CURRENT X') WRITE(*,'(" ")') RETURN END SUBROUTINE prt8 SUBROUTINE prt9(max, n, t, xopt, vm, fopt, nup, ndown, nrej, lnobds, nnew) REAL (dp), INTENT(IN) :: xopt(:), vm(:), t, fopt INTEGER, INTENT(IN) :: n, nup, ndown, nrej, lnobds, nnew LOGICAL, INTENT(IN) :: max Local variable INTEGER :: totmov totmov = nup + ndown + nrej WRITE(*,'(/," intermediate results before next temperature reduction",/)') WRITE(*,'(" CURRENT TEMPERATURE: ",G12.5)') t IF (max) THEN WRITE(*, '(" MAX FUNCTION VALUE SO FAR: ",G25.18)') fopt WRITE(*, '(" TOTAL MOVES: ",I8)') totmov WRITE(*, '(" UPHILL: ",I8)') nup WRITE(*, '(" ACCEPTED DOWNHILL: ",I8)') ndown WRITE(*, '(" REJECTED DOWNHILL: ",I8)') nrej WRITE(*, '(" OUT OF BOUNDS TRIALS: ",I8)') lnobds WRITE(*, '(" NEW MAXIMA THIS TEMPERATURE:",I8)') nnew ELSE WRITE(*, '(" MIN FUNCTION VALUE SO FAR: ",G25.18)') -fopt WRITE(*, '(" TOTAL MOVES: ",I8)') totmov WRITE(*, '(" DOWNHILL: ",I8)') nup WRITE(*, '(" ACCEPTED UPHILL: ",I8)') ndown WRITE(*, '(" REJECTED UPHILL: ",I8)') nrej WRITE(*, '(" TRIALS OUT OF BOUNDS: ",I8)') lnobds WRITE(*, '(" NEW MINIMA THIS TEMPERATURE:",I8)') nnew CALL prtvec(xopt, n, 'CURRENT OPTIMAL X') CALL prtvec(vm, n, 'STEP LENGTH (VM)') WRITE(*, '(" ")') RETURN END SUBROUTINE prt9 (110 of 114) 오전 10:21:36

111 SUBROUTINE prt10() WRITE(*, '(/, " SA ACHIEVED TERMINATION CRITERIA. IER = 0. ",/)') RETURN END SUBROUTINE prt10 SUBROUTINE prtvec(vector, ncols, name) This subroutine prints the double precision vector named VECTOR. Elements 1 thru NCOLS will be printed. NAME is a character variable that describes VECTOR. Note that if NAME is given in the call to PRTVEC, it must be enclosed in quotes. If there are more than 10 elements in VECTOR, 10 elements will be printed on each line. INTEGER, INTENT(IN) :: ncols REAL (dp), INTENT(IN) :: vector(ncols) CHARACTER (LEN=*), INTENT(IN) :: name INTEGER :: i, lines, ll WRITE(*,1001) NAME IF (ncols > 10) THEN lines = INT(ncols/10.) DO i = 1, lines ll = 10*(i - 1) WRITE(*,1000) vector(1+ll:10+ll) WRITE(*,1000) vector(11+ll:ncols) ELSE WRITE(*,1000) vector(1:ncols) 1000 FORMAT( 10(g12.5, ' ')) 1001 FORMAT(/, 25(' '), a) RETURN END SUBROUTINE prtvec END MODULE simulated_anneal PROGRAM simann USE simulated_anneal IMPLICIT NONE INTEGER, PARAMETER :: n = 2, neps = 4 REAL (dp) :: lb(n), ub(n), x(n), xopt(n), c(n), vm(n), t, eps, rt, fopt INTEGER :: ns, nt, nfcnev, ier, iseed1, iseed2, i, maxevl, iprint, & nacc, nobds LOGICAL :: max Set underflows to zero on IBM mainframes. CALL XUFLOW(0) (111 of 114) 오전 10:21:36

112 Set input parameters. max =.false. eps = 1.0D-6 rt =.5 iseed1 = 1 iseed2 = 2 ns = 20 nt = 5 maxevl = iprint = 1 DO i = 1, n lb(i) = -1.0D25 ub(i) = 1.0D25 c(i) = 2.0 Note start at local, but not global, optima of the Judge function. x(1) = x(2) = Set input values of the input/output parameters. t = 5.0 vm(1:n) = 1.0 WRITE(*,1000) n, max, t, rt, eps, ns, nt, neps, maxevl, iprint, iseed1, iseed2 CALL prtvec(x, n, 'STARTING VALUES') CALL prtvec(vm, n, 'INITIAL STEP LENGTH') CALL prtvec(lb, n, 'LOWER BOUND') CALL prtvec(ub, n, 'UPPER BOUND') CALL prtvec(c, n, 'C VECTOR') WRITE(*, '(/, " **** END OF DRIVER ROUTINE OUTPUT ****"/, & & " **** before CALL TO sa. ****")') CALL sa(n, x, max, rt, eps, ns, nt, neps, maxevl, lb, ub, c, iprint, iseed1, & iseed2, t, vm, xopt, fopt, nacc, nfcnev, nobds, ier) WRITE(*, '(/, " **** RESULTS AFTER SA **** ")') CALL prtvec(xopt, n, 'SOLUTION') CALL prtvec(vm, n, 'FINAL STEP LENGTH') WRITE(*,1001) fopt, nfcnev, nacc, nobds, t, ier 1000 FORMAT(/,' SIMULATED ANNEALING EXAMPLE',/,/, & ' NUMBER OF PARAMETERS: ',i3,' MAXIMIZATION: ',l5, /, & ' INITIAL TEMP: ', g8.2, ' RT: ',g8.2, ' EPS: ',g8.2, /, & ' NS: ',i3, ' NT: ',i2, ' NEPS: ',i2, /, & ' MAXEVL: ',i10, ' IPRINT: ',i1, ' ISEED1: ',i4, & ' ISEED2: ',i4) 1001 FORMAT(/,' OPTIMAL FUNCTION VALUE: ',g20.13 & /,' NUMBER OF FUNCTION EVALUATIONS: ',i10, & /,' NUMBER OF ACCEPTED EVALUATIONS: ',i10, & /,' NUMBER OF OUT OF BOUND EVALUATIONS: ',i10, & /,' FINAL TEMP: ', g20.13,' IER: ', i3) STOP END PROGRAM simann SUBROUTINE fcn(n, theta, h) This subroutine is from the example in Judge et al., The Theory and Practice of Econometrics, 2nd ed., pp There are two optima: F(.864,1.23) = (the global minumum) and F(2.35,-.319) = (112 of 114) 오전 10:21:36

113 IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN) :: theta(:) REAL (dp), INTENT(OUT) :: h Local variables INTEGER :: i REAL (dp) :: y(20), x2(20), x3(20) y(1) = y(2) = y(3) = y(4) = y(5) = y(6) = y(7) = y(8) = y(9) = y(10) = y(11) = y(12) = y(13) = y(14) = y(15) = y(16) = y(17) = y(18) = y(19) = y(20) = x2(1) =.286 x2(2) =.973 x2(3) =.384 x2(4) =.276 x2(5) =.973 x2(6) =.543 x2(7) =.957 x2(8) =.948 x2(9) =.543 x2(10) =.797 x2(11) =.936 x2(12) =.889 x2(13) =.006 x2(14) =.828 x2(15) =.399 x2(16) =.617 x2(17) =.939 x2(18) =.784 x2(19) =.072 x2(20) =.889 x3(1) =.645 x3(2) =.585 x3(3) =.310 x3(4) =.058 x3(5) =.455 x3(6) =.779 x3(7) =.259 x3(8) = (113 of 114) 오전 10:21:36

114 x3(9) =.028 x3(10) =.099 x3(11) =.142 x3(12) =.296 x3(13) =.175 x3(14) =.180 x3(15) =.842 x3(16) =.039 x3(17) =.103 x3(18) =.620 x3(19) =.158 x3(20) =.704 h = 0.0_dp DO i = 1, 20 h = (theta(1) + theta(n)*x2(i) + (theta(n)**2)*x3(i) - y(i))**2 + h RETURN END SUBROUTINE fcn (114 of 114) 오전 10:21:36