CS 404: Lab 3: Calling FORTRAN from CGetting Started1. Boot into Linux. To do this, click on “shutdown” and then select“restart.” The computer will begin restarting. You will eventuallyhave the option to select e ither Win2K or Linux (use up and downarrows to switch, press enter to select).2. Open a web browser and go to the course website1. Download RAD1Ditpack.tar.3. To unpack the tar-file, type tar xvf RAD1Ditpack.tar in a terminal win-dow. Cd into the directory “RAD1Ditpack.”Calling FORTRAN from C1. RAD1Ditpack contains the source code for simulating 1-D reaction-advection-diffusion problems with periodic boundary conditions. Thisis essentially the model-problem from CS403 with s ome slight modifi-cations2. As originally constructed, RAD1D solves for diffusion usingan implicit scheme. Computationally, this means that we have to solvea system of linear equations:A ∗ C = RHSfor each time step. The matrix A has several important properties, chiefamong them, that it is sparse (most entries are zeros) and symmetric(the lower triangular portion is a mirror image of the upper triangularportion, or mathematically, A(i, j) = A(j, i)). Because A is sparse,we can save a tremendous amount of memory by only tracking thelocations where A 6= 0 rather than storing all of the zeros. If n is thenumber of rows in A, then storing A as a full 2D array would take8n2bytes. Since there are only 3n nonzero elements, we can reducethe memory requirements to 24n bytes plus an additional 16n bytesfor the integer arrays needed to remember where A is nonzero. This iswhat the original RAD1D did. Because A is symmetric, we can further1www.cs.cornell.edu/Course/cs404/2002sp2Check out www.cs.cornell.edu/Courses/cs403/Lecture08/note.html for more info onRAD1D1reduce the storage requirements by only storing the upper triangularportion of A. I’ve modified main.c and the routine BuildA in linalg.cto only store the nonzeros in the upper triangular portion. Take a lookif you’re interested.2. Unfortunately, the routine to solve Ac = b in the original version (pcgmin linalg.c) requires every nonzero entry in A to be represented. Per-haps more unfortunately, to use the routines in LAPACK with thisproblem, we would have to represent the entire upper triangular por-tion of A–both zeros and nonzeros. In order to take advantage ofA’s sparsity and symmetry, we need a new solver. Going through theGAMS tree, I found the package ITPACK which provides a variety ofroutines for solving sparse linear systems. The only problem: ITPACKis in FORTRAN, and RAD1D is in C, so we need to figure out how tocall ITPACK from our C code. After reading the ITPACK documen-tation (a PDF version is on the course web site), I figured out that weneed to call two functions: DFAULT which places the default parame-ters into two arrays IPARM and RPARM (integer and real parameters,respectively), and JCG which solves our system of equations using the“conjugate gradient method with a Jacobi preconditioner.” I have re-placed the call to the original pcgm routine (in SolveA in linalg.c) withcalls to ITPACK, and I have also allocated the workspace arrays re-quired by JCG (WKSP and IWKSP) and set the length of these arrays(NW). Your job is to get the calls to DFAULT and JCG to work. How-ever, rather than talk you through it step-by-step as I usually do, I willonly give you a checklist of things to do:2Key Description1. Call-by-value Make sure all variables are called by reference (areeither arrays or pointers) rather than by value.Remember, the reference to a C variable is ob-tained with the & operator: foo–value of variablefoo, &foo–memory address (reference) of variablefoo2. Case Most FORTRAN compilers (including g77 makeall variables and function calls lower case; if theydon’t you can usually force everything lower casewith a compiler option (usually -f). For every-thing to link correctly, you must call FORTRANroutines using a lower case name.3. Underscore Most FORTRAN compilers append an underscore(“ ”) to each subroutine name. For everything tolink correctly, you may need to add the underscorethe subroutine call in the C program.4. Prototype Strict ANSI C requires you to provide a proto-type for every subroutine. The prototype givesthe return type of the subroutine (void, if no val-ues are explicitly returned), subroutine name, andthe types of each variable. You need to provideprototypes for both DFAULT and JCG. There’sa special “FORTRAN prototype” section at thetop of linalg.c3. When you think you’ve made the changes in linalg.c, you’re ready tobuild the executable. I’ve modified the Makefile to compile the IT-PACK routines in dsrc2c.f to dsrc2c.o. Typing “make” should createall of the objects and attempt to link them using f77. Why do we linkwith f77 rather than gcc? If you made the right changes in SolveA,rad1d should build correctly. You can run the problem by typing rad1dbasic.cmnd (basic.cmnd is a “command” file that specifies the parame-ters for the simulation–see the online description for more