Initial commit with code and start of report

9 files changed, 570 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..484ab7e
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+results/*
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..8c48932
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,18 @@
+FC 		= mpiifx
+#FFLAGS 	= -o $(TARGET) -qmkl=sequential
+FFLAGS 	= -o $(TARGET) -qmkl
+TARGET 	= test.x
+NWORK 	= 64
+
+all: clean $(TARGET) run
+
+$(TARGET): main.f90 worker.f90 manager.f90
+	$(FC) $(FFLAGS) main.f90 worker.f90 manager.f90
+
+clean:
+	rm -f *.x *.mod eigs
+
+run:
+	mpirun ./$(TARGET) -np $(NWORK)
+
+
diff --git a/main.f90 b/main.f90
new file mode 100644
index 0000000..e62d47f
--- /dev/null
+++ b/main.f90
@@ -0,0 +1,117 @@
+include 'mkl_vsl.f90'                     ! For quasi-random number generation with the MKL VSL library
+module globals
+USE MKL_VSL_TYPE                          ! Routines for quasi-random 
+USE MKL_VSL                               ! number generation
+implicit none
+include 'mpif.h'                          ! Fortran MPI header file
+integer :: n                              ! Size of the matrices
+integer, parameter :: EXIT_TAG=1,DEFAULT_TAG=0,void=0 ! Used to signal workers to return
+integer :: proc_num,num_procs             ! Process number and number of processes
+end module globals
+
+module auxiliary
+use globals
+implicit none
+TYPE (VSL_STREAM_STATE) :: stream         ! Identifier for the pseudo-random number stream
+contains
+subroutine get_random_seed(seed)
+integer :: un,seed,istat
+
+! Read seed integer from /dev/urandom, a file filled with pseudo-random numbers generated from the state of the computer
+open(newunit=un, file="/dev/urandom", access="stream", &
+                 form="unformatted", action="read", status="old", iostat=istat)
+   read(un) seed
+close(un)
+
+end subroutine get_random_seed
+
+subroutine open_random_stream(seed)
+! Setup of stream of pseudo-random numbers using MKL library.
+integer :: ierr,seed
+
+! Use a Mersenne Twister pseudorandom number generator; the PRN sequence is identified by the variable "stream"
+ierr=vslnewstream(stream,VSL_BRNG_MT19937,seed)
+if(ierr.ne.0) then
+   print *,'Trouble with the MKL RNG, returns flag ',ierr
+end if
+
+end subroutine open_random_stream
+
+subroutine close_random_stream
+! Close pseudo-random number stream
+integer :: rnd_ierr
+
+rnd_ierr=vsldeletestream( stream )
+
+end subroutine close_random_stream
+end module auxiliary
+
+program main
+use globals
+use auxiliary
+implicit none
+logical :: ok
+double precision :: wtime
+
+call start_MPI(ok)
+if(proc_num.eq.0) wtime=MPI_wtime()
+
+if(ok) then
+   if(proc_num.eq.0) then
+      call manager
+   else
+      call worker
+   end if
+end if
+if(proc_num.eq.0) then
+   wtime=MPI_wtime()-wtime
+   print *,'wtime=',wtime,'(s)'
+end if
+call stop_MPI
+
+end program main
+
+subroutine start_MPI(ok)
+use globals
+implicit none
+! MPI initialization
+logical :: my_ok=.true.,ok
+integer :: ierr
+call mpi_init(ierr)
+
+if(ierr.ne.0) then
+   print *,'MPI_init failed!'
+   my_ok=.false.
+end if
+
+if(my_ok) then
+   call mpi_comm_size(MPI_COMM_WORLD,num_procs,ierr)
+   if(ierr.ne.0) then
+      print *,'MPI_comm_size failed!'
+      my_ok=.false.
+   end if
+   if(my_ok) then
+      call mpi_comm_rank(MPI_COMM_WORLD,proc_num,ierr)
+      if(ierr.ne.0) then
+         print *,'MPI_comm_rank failed!'
+         my_ok=.false.
+      end if
+   end if
+end if
+! Check if everyone is ok.
+call mpi_allreduce(my_ok,ok,1,MPI_LOGICAL,MPI_LAND,MPI_COMM_WORLD,ierr)
+end subroutine start_MPI
+subroutine stop_MPI
+use globals
+implicit none
+logical :: init
+integer :: ierr
+
+! Wait until everybody is done. One process "finalize"ing while others are still working can cause ugly crashes.
+call mpi_barrier(MPI_COMM_WORLD,ierr)
+! Check if MPI has been initialized.
+call mpi_initialized(init,ierr)
+! If it is, call finalize.
+if(init) call mpi_finalize(ierr)
+
+end subroutine stop_MPI
diff --git a/manager.f90 b/manager.f90
new file mode 100644
index 0000000..de85dde
--- /dev/null
+++ b/manager.f90
@@ -0,0 +1,79 @@
+subroutine manager
+use globals
+use auxiliary
+implicit none
+logical :: start(1:num_procs-1)
+integer :: seed,ierr,recvd=0,worker,exit_tags_sent,tag,ndat,p,failed
+integer, dimension(mpi_status_size) :: status
+double precision :: buf
+double precision, allocatable, dimension(:) :: eigs
+
+! Set matrix size and number of eigenvalues
+call init(ndat)
+! Broadcast matrix size
+call mpi_bcast(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
+! Allocate array for received eigenvalues
+allocate(eigs(ndat))
+
+start=.true.
+exit_tags_sent=0
+failed=0
+open(22,file='eigs')
+write(22,'(i0,i0)') n, ndat
+
+! Receive data until at least ndat have been received
+do while(.true.)
+   ! For each worker check if a seed is needed
+   do p=1,num_procs-1
+      if(start(p)) then
+         call get_random_seed(seed)
+         if(recvd.lt.ndat-num_procs+2) then
+            tag=DEFAULT_TAG
+         else
+            tag=EXIT_TAG
+            exit_tags_sent=exit_tags_sent+1
+         end if
+         call mpi_send(seed,1,MPI_INTEGER,p,tag,MPI_COMM_WORLD,ierr)
+!         print *,'Manager sent seed ',seed,' to worker ',p
+         start(p)=.false.
+      end if
+   end do
+   if(exit_tags_sent.eq.num_procs-1) then
+      print *,'Manager exiting...'
+      exit
+   end if
+   
+   call mpi_recv(buf,1,MPI_DOUBLE_PRECISION,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,status,ierr)
+   worker = status(MPI_SOURCE)
+   tag = status(MPI_TAG)
+   if(tag.eq.0) then
+      recvd=recvd+1
+      eigs(recvd)=buf
+      write(22,'(i4,e14.7)') worker,buf
+   else
+      failed=failed+1
+      print '(i5,a7,i5,a8)',recvd,' received, ',failed,' failed'
+   end if
+   start(worker)=.true.
+end do
+
+close(22)
+
+deallocate(eigs)
+return
+end subroutine manager
+
+subroutine init(ndat)
+use globals
+implicit none
+integer :: ndat
+
+! Matrix dimension and sample size.. hard-coded now, but can be read from command line or file...
+!read(*,*) n
+!read(*,*) ndat
+
+n=4096
+ndat=2048
+
+return
+end subroutine init
diff --git a/report/.gitignore b/report/.gitignore
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/report/.gitignore
diff --git a/report/refs.bib b/report/refs.bib
new file mode 100644
index 0000000..cf0a85f
--- /dev/null
+++ b/report/refs.bib
@@ -0,0 +1,82 @@
+@book{novak2022,
+    year = {2022},
+    author = {Kyle A. Novak},
+    isbn = {9798985421804},
+    title = {Numerical Methods for Scientific Computing},
+    note = {Chapter 4},
+    month = {March},
+    publisher = {Equal Share Press},
+    address = {Sharon, Vermont},
+}
+
+@book{watkins2002,
+    title = {Fundamenals of Matrix Computations (2nd Edition)},
+    month = {January},
+    address = {New York},
+    year = {2002},
+    author = {David S. Watkins},
+    isbn = {0471213942},
+    note = {Chapter 5},
+    publisher = {John Wiley and Sons},
+}
+
+@article{golub2009The_QR,
+    pages = {467-485},
+    issn = {1464-3642},
+    number = {3},
+    doi = {10.1093/imanum/drp012},
+    publisher = {Oxford University Press (OUP)},
+    volume = {29},
+    title = {{The QR Algorithm: 50 Years Later its Genesis by John Francis and Vera Kublanovskaya and Subsequent Developments}},
+    journal = {IMA Journal of Numerical Analysis},
+    month = {June},
+    author = {Golub, G. and Uhlig, F.},
+    note = {},
+    url = {http://dx.doi.org/10.1093/imanum/drp012},
+    year = {2009},
+}
+
+@article{wigner1955Charac,
+    number = {3},
+    publisher = {JSTOR},
+    pages = {548},
+    url = {http://dx.doi.org/10.2307/1970079},
+    title = {{Characteristic Vectors of Bordered Matrices With Infinite Dimensions}},
+    doi = {10.2307/1970079},
+    journal = {The Annals of Mathematics},
+    issn = {0003-486X},
+    year = {1955},
+    volume = {62},
+    author = {Wigner, Eugene P.},
+}
+
+@book{anderson2009An_Int,
+    isbn = {9781107471580},
+    author = {Anderson, Greg W. and Guionnet, Alice and Zeitouni, Ofer},
+    month = {February},
+    publisher = {Cambridge University Press},
+    doi = {10.1017/cbo9780511801334},
+    url = {http://dx.doi.org/10.1017/CBO9780511801334},
+    note = {Chapters 1 and 2},
+    title = {{An Introduction to Random Matrices}},
+    year = {2009},
+}
+
+@misc{marco2013,
+    author = {Marco},
+    title = {{Approximation for the Tracy-Widom Laws}},
+    year = {2013},
+    howpublished = {Published to the MATLAB Central File Exchange},
+    note = {Retrieved April 25th, 2026. \href{https://www.mathworks.com/matlabcentral/fileexchange/44711-approximation-for-the-tracy-widom-laws}{Available online.}},
+}
+
+@article{chiani2012Distri,
+    doi = {10.48550/ARXIV.1209.3394},
+    url = {https://arxiv.org/abs/1209.3394},
+    author = {Chiani, Marco},
+    keywords = {Information Theory (cs.IT), Statistics Theory (math.ST), FOS: Computer and information sciences, FOS: Computer and information sciences, FOS: Mathematics, FOS: Mathematics},
+    title = {{Distribution of the Largest Eigenvalue for Real Wishart and Gaussian Random Matrices and a Simple Approximation for the Tracy-Widom Distribution}},
+    publisher = {arXiv},
+    copyright = {arXiv.org perpetual, non-exclusive license},
+    year = {2012},
+}
\ No newline at end of file
diff --git a/report/report.tex b/report/report.tex
new file mode 100644
index 0000000..7b66883
--- /dev/null
+++ b/report/report.tex
@@ -0,0 +1,101 @@
+\documentclass{article}
+
+\usepackage[margin=1in]{geometry}
+\usepackage{parskip}
+\usepackage{float}
+
+\usepackage{xcolor}
+\usepackage{graphicx}
+\graphicspath{{./figures/}}
+\usepackage{amsmath}
+\usepackage{booktabs}
+\usepackage{longtable}
+
+\usepackage{hyperref}
+\hypersetup{colorlinks=true,
+            linkcolor=black,
+            urlcolor=blue,
+            citecolor=black,
+            pdftitle={Reconstructing the Tracy-Widom Distribution},
+            pdfpagemode=FullScreen}
+
+\usepackage{listings}
+\usepackage{algpseudocode}
+
+\title{Reconstructing the Tracy-Widom Distribution\vspace{-0.75em}}
+\author{Connor Moore, 100826701. \today{}}
+\date{}
+
+\begin{document}
+\maketitle
+
+\section{Introduction}
+Random numbers form the basis of various studies in mathematics and the sciences. One strong example of this is \emph{Random Matrix Theory}, or the study of the construction and properties of random matrices. One of the first applications of random matrix theory was proposed by Eugene Wigner (1902-1995) in his application to model the nuclei and spectrum of various heavy atom \cite{anderson2009An_Int, wigner1955Charac}. It has since found different applications including work in computational mechanics, control theory, systems engineering, and others.
+
+Random matrix theory borrows some distributions from statistics to construct matrices. This includes the Gaussian (normal) distribution, which is used to make various matrices depending on which number system is used. Wigner's work included Orthogonal (real symmetric), Unitary (Complex Hermitian), and Symplectic (Quaternion Hermitian) matrices and their eigenvalues. 
+
+This work considers Gaussian orthogonal ensembles (GOEs) specifically and attempts to reconstruct the Tracy-Widom distribution directly by computing eigenvalues. An overview of relevant theory, the necessary methods, and implementation details is provided before results are presented. The distribution is examined as a function of matrix dimension and sample size, and statistical testing is performed to support the notion of a Tracy-Widom being obtained.
+
+\section{Relevant Theory}
+
+\subsection{Random Ensembles}
+The GOE-based random matrix has some properties that go beyond simply being normally distributed. Namely, the matrix is square, symmetric, and the standard deviation $\sigma$ is unity except for elements on the diagonal where $\sigma=2$. The mean $\mu$ is considered as zero for all elements. Consider the normal distribution with notation $N(\mu, \sigma$). Then, an example is provided below for a $5\times5$ matrix:
+\begin{equation*}
+    \mathbf{A} = 
+    \begin{pmatrix}
+        N(0,2) & N(0,1) & N(0,1) & N(0,1) & N(0,1) \\
+        N(0,1) & N(0,2) & N(0,1) & N(0,1) & N(0,1) \\
+        N(0,1) & N(0,1) & N(0,2) & N(0,1) & N(0,1) \\
+        N(0,1) & N(0,1) & N(0,1) & N(0,2) & N(0,1) \\
+        N(0,1) & N(0,1) & N(0,1) & N(0,1) & N(0,2) \\
+    \end{pmatrix}
+\end{equation*}
+The diagonal standard deviation of 2 is what makes this ensemble orthogonal. Note that the representation above does not explicitly state the symmetry, but $\mathbf{A}(i,j) = \mathbf{A}(j,i)\, \forall(i,j)$.
+
+\subsection{The Tracy-Widom Distributions}
+Various distributions can be 
+
+\subsection{Statistical Differentiation}
+The Tracy-Widom and normal distributions are hard to differentiate based on shape alone. Because of this, it is necessary to consider quantifying the differences in some meaningful way. One way to do this is to consider the important differences between the two distributions. For example, while the normal distribution is symmetric, the Tray-Widom is not. The symmetry of a distribution is something that can be quantified
+
+\section{Numerical Methods}
+\subsection{The Power Method}
+One relatively cheap way to find the maximum eigenvalue is through the so-called power method. The power method is an iterative approach to finding the maximum eigenvalue of a matrix that exploits the spectral gap ($\lambda_1 - \lambda_2$). Consider an initial guess in the form of a vector of length $N$ given as $\mathbf{x_0}$. Then, perform the multiplication:
+\begin{equation}
+    \mathbf{x}_{k+1} = \mathbf{A}\mathbf{x}_k
+\end{equation}
+After multiplying, then normalize the vector such that $\left| \mathbf{x}_{k+1} \right|$ is unity:
+\begin{equation}
+    \mathbf{x}_{k+1} = \frac{\mathbf{x}_{k+1}}{\left| \mathbf{x}_{k+1} \right|}
+\end{equation}
+Or, more compactly written as one step:
+\begin{equation}
+    \mathbf{x}_{k+1} = \frac{\mathbf{A}\mathbf{x}_{k}}{\left|\mathbf{A}\mathbf{x}_{k}\right|}
+\end{equation}
+Repeating this will converge to the dominant eigenvalue $\lambda_1$ so long as $\mathbf{x}_0$ is not orthogonal to the dominant eigenvector \cite{novak2022}. The convergence of the method is $\mathcal{O}\left( \left| \lambda_1 / \lambda_2 \right|^k \right)$ and is largely dependent on the spectral gap. The main advantage of this method is that computationally it consists of only a matrix multiplication ($\mathcal{O}(n^2)$) and normalization ($\mathcal{O}(n)$) which is rather inexpensive compared with other methods per iteration \cite{watkins2002}. Worth noting, however, is that in cases where $\mathbf{A}$ has a small spectral gap it is not guaranteed to converge in an acceptable number of iterations. This poses an issue as it would bias the calculations for reconstructing the Tracy-Widom by not including random ensembles specifically with lower spectral gaps. To ratify this, more expensive methods exist which guarantee a solution that are used as a fallback.
+\subsection{Schur Decomposition}
+One method that guarantees a solution $\forall \mathbf{A}$ is Schur Decomposition. Named after the late Russian mathematician Issai Schur (1875-1941), Schur Decomposition involves decomposing the matrix $\mathbf{A}$ into a combination of some unitary matrix $\mathbf{U}$ and upper-triangular matrix $\mathbf{T}$ \cite{watkins2002}:
+\begin{equation}
+    \mathbf{A} = \mathbf{U} \mathbf{T} \mathbf{U}^\mathrm{T}
+\end{equation}
+In this case, the diagonal entries in $\mathbf{T}$ are then the eigenvalues of the original matrix $\mathbf{A}$. The proof in Schur's theorem is non-constructive and does not provide a solution without knowing the eigenvectors of $\mathbf{A}$, but it can be still be used in combination with the so-called QR algorithm developed independently by John Francis (1934-present) and Vera Kublanovskaya (1920-2012) starting from 1959 \cite{golub2009The_QR}. Essentially, the QR algorithm is the computational scheme that provides the matrices $\mathbf{U}$ and $\mathbf{T}$ from the input $\mathbf{A}$. It does so with a complexity of $\mathcal{O}(n^3)$, which for larger values of $N$ deviates significantly from the $\mathcal{O}(n^2)$ of the power method. For this reason it is used only when the power method fails.
+
+\section{Implementation and MPI Parallelization}
+
+\subsection{Worker-Manager Structure}
+
+\subsection{The Manager Process}
+
+\subsection{The Calculation (Worker) Process}
+
+
+\section{Results}
+
+%\section{Conclusions}
+
+\bibliographystyle{ieeetr}
+\bibliography{refs.bib}
+
+\end{document}
+    
+
diff --git a/run_sweeps.sh b/run_sweeps.sh
new file mode 100755
index 0000000..ebb229b
--- /dev/null
+++ b/run_sweeps.sh
@@ -0,0 +1,43 @@
+#!/bin/bash
+
+
+#N_MATRIX="1024 2048 4096 8192 16384"
+N_MATRIX="1024 2048 4096"
+N_DATAPTS="64 128 256 512 1024 2048"
+
+GREEN='\033[0;32m'
+RESET='\033[0m'
+
+mkdir -p results/logs
+touch results/wtable
+
+
+for N in $N_MATRIX; do
+	# Change the manager file to use the correct matrix size
+	echo -e "${GREEN}Updating matrix size to $N...${RESET}"
+	sed "s/n=[0-9]\+/n=$N/" manager.f90 -i
+	for L in $N_DATAPTS; do
+		# Change the worker file to use the correct number of datapts
+		echo -e "Setting data size to $L!"
+		sed "s/ndat=[0-9]\+/ndat=$L/" manager.f90 -i
+
+		# Print to double check
+		cat manager.f90 | grep "n="
+		cat manager.f90 | grep "ndat="
+
+		# And compile the new binary + run
+		LOGFILE=results/logs/$N\_$L.log
+		make > $LOGFILE
+
+		WTIME=$(cat $LOGFILE | grep "wtime=" | cut -d "=" -f 2 | cut -d "(" -f 1)
+		echo -e "Wall time was $WTIME\n"
+
+		# Append to wall time file
+		echo $N $L $WTIME >> results/wtable
+
+
+		# Move the results to the proper file
+		mv eigs results/eigs_$N\_$L.out
+
+	done
+done
diff --git a/worker.f90 b/worker.f90
new file mode 100644
index 0000000..08fdfa2
--- /dev/null
+++ b/worker.f90
@@ -0,0 +1,129 @@
+subroutine worker
+use globals
+use auxiliary
+implicit none
+logical :: conv
+integer :: seed,ierr,recvd_tag,tag
+integer, dimension(mpi_status_size) :: status
+double precision :: eig
+double precision, allocatable, dimension(:,:) :: A
+
+! Receive matrix size
+call mpi_bcast(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
+allocate(A(n,n))
+
+do while(.true.)
+   call mpi_recv(seed,1,MPI_INTEGER,0,MPI_ANY_TAG,MPI_COMM_WORLD,STATUS,ierr)
+   recvd_tag = status (MPI_TAG)
+   if(recvd_tag.eq.EXIT_TAG)  then            ! The exit tag was received, wind down and exit.
+      print *,'Worker ',proc_num,' exiting...'
+      exit
+   end if
+   ! The seed was received and more data are requested, so compute an eigenvalue:
+   call open_random_stream(seed)
+   call make_random_matrix(A,seed)
+   call compute_eig(A,eig,conv)
+   ! Send result
+   if(conv) then
+      tag=0
+   else
+      tag=1
+   end if
+   call mpi_send(eig,1,MPI_DOUBLE_PRECISION,0,tag,MPI_COMM_WORLD,ierr)
+   call close_random_stream()
+end do
+
+deallocate(A)
+return
+end subroutine worker
+
+subroutine make_random_matrix(A,seed)
+use globals
+use auxiliary
+implicit none
+integer :: seed,ierr,i,j
+double precision :: A(n,n)
+
+do i=1,n
+   ierr=vdRNGGaussian(VSL_RNG_METHOD_GAUSSIAN_ICDF,stream,1,A(i,i),0d0,2d0)
+   do j=1,i-1
+      ierr=vdRNGGaussian(VSL_RNG_METHOD_GAUSSIAN_ICDF,stream,1,A(i,j),0d0,1d0)
+      A(j,i)=A(i,j)
+   end do
+end do
+
+return
+end subroutine make_random_matrix
+
+subroutine compute_eig(A,eig,conv)
+use globals
+use auxiliary
+implicit none
+logical :: conv
+integer :: ierr,M,lwork,info,i
+double precision :: A(n,n),eig,eps,mu,wr(n),wi(n),B(n,n),vl(1,n),vr(1,n),eig2
+double precision, allocatable, dimension(:) :: work
+
+lwork=4*n
+allocate(work(lwork)) ! Bit of memory needed by LAPACK
+
+eps=1d-5              ! Relative convergence criterion for power iteration (bad practice to hard-code this!)
+M=int(n/2)            ! Maximal number of iterations
+mu=0d0
+call power(A,mu,eps,M,eig,conv)
+! Check if the eigenvalue is positive (they are symmetrically distributed about 0). If not, re-compute with shift:
+if(conv) then
+   if(eig.lt.0) then
+      mu=-eig
+      call power(A,mu,eps,M,eig,conv)
+   end if
+end if
+
+! If power iteration fails to converge, compute the whole spectrum with Shur decomposition:
+if(.not.(conv)) then
+  call DGEEV('N','N',n,A,n,wr,wi,vl,1,vr,1,work,lwork,info)
+  if(info.eq.0) then
+     eig=wr(maxloc(wr,1)) ! Select the largest eigenvalue
+     conv=.true.
+  else
+     print *,'error in DGEEV ',info
+  end if
+end if
+
+deallocate(work)
+
+return
+end subroutine compute_eig
+
+subroutine power(A,mu,eps,M,eig,conv)
+use globals
+use auxiliary
+implicit none
+logical :: conv
+integer :: i,M,ierr,j
+double precision :: A(n,n),v(n),y(n),eig,eps,mu,err,nrm,ndp,mem
+
+conv=.false.
+
+! Make random initial vector
+ierr=vdrnggaussian(VSL_RNG_METHOD_GAUSSIAN_ICDF,stream,n,v,0d0,1d0)
+v=v/sqrt(sum(v**2))
+mem=-100d0
+
+do i=1,M
+   call dgemv('N',n,n,1d0,A,n,v,1,0d0,y,1)
+   y=y+mu*v
+   nrm=sqrt(sum(y**2))
+   eig=dot_product(v,y) ! Rayleigh quotient approximates the eigenvalue
+   eig=eig-mu
+   err=abs(mem-eig)/abs(mem)
+   if(err.lt.eps) then
+      conv=.true.
+      exit
+   end if
+   v=y/nrm
+   mem=eig
+end do
+
+return
+end subroutine power