|
|
@@ -0,0 +1,1138 @@
|
|
|
+/******************************************************************************
|
|
|
+ S L A L O M
|
|
|
+
|
|
|
+ Scalable Language-independent Ames Laboratory One-minute Measurement
|
|
|
+
|
|
|
+ The following program is the first benchmark based on fixed time rather
|
|
|
+ than fixed problem comparison. Not only is fixed time more representative
|
|
|
+ of the way people use computers, it also greatly increases the scope and
|
|
|
+ longevity of the benchmark. SLALOM is very scalable, and can be used to
|
|
|
+ compare computers as slow as 126 floating-point operations per second
|
|
|
+ (FLOPS) to computers running a trillion times faster. The scalability can
|
|
|
+ be used to compare single processors to massively parallel collections
|
|
|
+ of processors, and to study the space of problem size vs. ensemble size
|
|
|
+ in fine detail. It resembles the LINPACK benchmark since it involves
|
|
|
+ factoring and backsolving a (nearly) dense matrix, but incorporates a
|
|
|
+ number of improvements to that benchmark that we hope will make SLALOM
|
|
|
+ a better reflection of general system performance.
|
|
|
+
|
|
|
+ The SLALOM benchmark solves a complete, real problem (optical radiosity
|
|
|
+ on the interior of a box), not a contrived kernel or a synthetic mixture of
|
|
|
+ sample operations. SLALOM is unusual since it times input, problem setup,
|
|
|
+ solution, and output, not just the solution. For slower computers, the
|
|
|
+ problem setup will take the majority of the time; it grows as the square of
|
|
|
+ the problem size. The solver grows as the cube of the problem size, and
|
|
|
+ dominates the time for large values of n.
|
|
|
+
|
|
|
+ While the following is C, you are free to translate it into any
|
|
|
+ language you like, including assembly language specific to one computer.
|
|
|
+ You may use compiler directives, hand-tuned library calls, loop unrolling,
|
|
|
+ and even change the algorithm, if you can provide a convincing argument
|
|
|
+ that the program still works for the full range of possible inputs. For
|
|
|
+ example, if you replace the direct solver with an iterative one, you must
|
|
|
+ make sure your method is correct even when the geometry is quite eccentric
|
|
|
+ and the box faces are highly reflective. (rho = .999)
|
|
|
+
|
|
|
+ The main() driver should be used with the value of 60 seconds for the
|
|
|
+ SLALOM benchmark. The work done for a particular problem size is figured
|
|
|
+ after timing has ceased, so there is no overhead for work assessment. The
|
|
|
+ residual check ||Ax - b|| is also done after timing has ceased. Two
|
|
|
+ computers may be compared either by their problem size n, or by their MFLOPS
|
|
|
+ rate, never by the ratio of execution times. Times will always be near one
|
|
|
+ minute in SLALOM. We have used the following weights for floating-point
|
|
|
+ operation counting, based on the weights used by Lawrence Livermore National
|
|
|
+ Laboratory:
|
|
|
+
|
|
|
+ OPERATION WEIGHT
|
|
|
+ a=b, a=(constant) 0
|
|
|
+ a<0, a<=0, a==0, a!=0, a>0, a>=0 0
|
|
|
+ -a, fabs(a), fsgn(a, b) 0
|
|
|
+ a+b, a-b, a*b, a^2 1
|
|
|
+ a<b, a<=b, a==b, a!=b, a>b, a>=b 1
|
|
|
+ (int) a, (double)b 1
|
|
|
+ 1/a, -1/a 3
|
|
|
+ a/b 4
|
|
|
+ sqrt(a) 4
|
|
|
+ Format to or from ASCII string 6
|
|
|
+ sin(a), cos(a), tan(a), log(a), atan(a), exp(a) 8
|
|
|
+
|
|
|
+ We invite you to share with us the results of any measurements that you
|
|
|
+ make with SLALOM. We do NOT accept anonymous data; machine timings will be
|
|
|
+ referenced and dated.
|
|
|
+
|
|
|
+ The least you need to do to adapt SLALOM to your computer is:
|
|
|
+
|
|
|
+ 1. In the "Measure" routine, set NMAX to a value large enough to keep
|
|
|
+ the computer working for a minute. Vary it slightly if it helps
|
|
|
+ (for reasons of cache size, interleaving, etc.)
|
|
|
+
|
|
|
+ 2. Replace the timer call in "When" with the most accurate wall-clock
|
|
|
+ timer at your disposal. If only CPU time is available, try to run
|
|
|
+ the job standalone or at high priority, since we are ultimately
|
|
|
+ interested in the top of the statistical range of performance.
|
|
|
+
|
|
|
+ 3. Edit in the information specific to your test in the "What"
|
|
|
+ routine, so that final output will be automatically annotated.
|
|
|
+
|
|
|
+ 4. Compile, link, and run the program, interacting to select values
|
|
|
+ of n that bracket a time of one minute. Once everything is
|
|
|
+ running, run it as a batch job so as to record the session.
|
|
|
+
|
|
|
+ Examples of ways you may optimize performance:
|
|
|
+
|
|
|
+ 1. Unroll the loops in SetUp1 and SetUp2; it is possible to
|
|
|
+ vectorize both SetUp1 and SetUp2 at the cost of some extra
|
|
|
+ operations, program complexity, and storage.
|
|
|
+
|
|
|
+ 2. Replace the innermost loops of Solver with calls to well-tuned
|
|
|
+ libraries of linear algebra routines, such as DDOT from the
|
|
|
+ Basic Linear Algebra Subroutines (level 1 BLAS). Better still,
|
|
|
+ use a tuned library routine for all of Solver; the sparsity
|
|
|
+ exploited in Solver is only a few percent, so you will usually
|
|
|
+ gain more than you lose by applying a dense symmetric solver.
|
|
|
+
|
|
|
+ 3. Parallelize the SetUp and Solver routines; all are highly
|
|
|
+ parallel. Each element of the matrix can be constructed
|
|
|
+ independently, once each processor knows the geometry and part of
|
|
|
+ the partitioning into regions. A substantial body of literature
|
|
|
+ now exists for performing the types of operations in Solver in
|
|
|
+ parallel.
|
|
|
+
|
|
|
+ 4. Overlap computation with output. Once the Region routine is done,
|
|
|
+ the first part of the output file (patch geometry) can be written
|
|
|
+ while the radiosities are being calculated.
|
|
|
+
|
|
|
+ Examples of what you may NOT do:
|
|
|
+
|
|
|
+ 1. The tuning must not be made specific to the particular input
|
|
|
+ provided. For example, you may not eliminate IF tests simply
|
|
|
+ because they always come out the same way for this input; you
|
|
|
+ may not use precomputed answers or table look-up unless those
|
|
|
+ answers and tables cover the full range of possible inputs; and
|
|
|
+ you may not exploit symmetry for even values of the problem size.
|
|
|
+
|
|
|
+ 2. You may not disable the self-consistency tests in SetUp3 and
|
|
|
+ Verify, nor alter their tolerance constants.
|
|
|
+
|
|
|
+ 3. You may not change the input or output files to unformatted
|
|
|
+ binary or other format that would render them difficult to create
|
|
|
+ or read for humans.
|
|
|
+
|
|
|
+ 4. You may not eliminate the reading of the "geom" file by putting
|
|
|
+ its data directly into the compiled program.
|
|
|
+
|
|
|
+ 5. You may not change any of the work assessments in Meter. If you
|
|
|
+ use more floating-point operations than indicated, you must still
|
|
|
+ use the assessments provided. If you find a way to use fewer
|
|
|
+ operations and still get the job done for arbitrary input
|
|
|
+ parameters, please tell us!
|
|
|
+
|
|
|
+ -John Gustafson, Diane Rover, Michael Carter,
|
|
|
+ and Stephen Elbert
|
|
|
+ Ames Laboratory, Ames, Iowa 50011
|
|
|
+******************************************************************************/
|
|
|
+
|
|
|
+/*****************************************************************************/
|
|
|
+/* The following program finds a value n such that a problem of size n */
|
|
|
+/* takes just under "goal" seconds to execute. */
|
|
|
+/* */
|
|
|
+/* John Gustafson, Diane Rover, Michael Carter, and Stephen Elbert */
|
|
|
+/* Ames Laboratory, 3/18/90 */
|
|
|
+/* */
|
|
|
+/* Calls: Meter Measures execution time for some application. */
|
|
|
+/* What Prints work-timing statistics and system information. */
|
|
|
+/*****************************************************************************/
|
|
|
+
|
|
|
+#include <stdio.h>
|
|
|
+#include <math.h>
|
|
|
+#include <sys/time.h>
|
|
|
+
|
|
|
+/* NMAX = Largest npatch for your computer; adjust as needed. */
|
|
|
+#define NMAX 2048
|
|
|
+#define EPS (0.5e-8)
|
|
|
+#define FALSE (1==0)
|
|
|
+#define TRUE (!FALSE)
|
|
|
+#define MAX(a,b) (((a) > (b)) ? (a) : (b))
|
|
|
+
|
|
|
+/* Global variables and function return types: */
|
|
|
+double goal, /* User input, fixed-time benchmark goal, in seconds. */
|
|
|
+ timing, /* Elapsed time returned by Meter routine, in seconds.*/
|
|
|
+ work, /* In this case, number of FLOPs performed. */
|
|
|
+ When(), /* Wall clock in seconds. */
|
|
|
+ Ddot(); /* Double dot product. */
|
|
|
+int mean, /* Avg between upper and lower bounds for bisection */
|
|
|
+ /* method. */
|
|
|
+ n, /* The problem size. */
|
|
|
+ nupper, /* Upper bound on problem size, used in iterating */
|
|
|
+ /* toward goal. */
|
|
|
+ Meter(), /* Driver for following benchmark functions. */
|
|
|
+ Reader (), /* Reads problem description from 'geom' file. */
|
|
|
+ Region (), /* Subdivides box faces into patches. */
|
|
|
+ SetUp3 (), /* Set up matrix to solve. */
|
|
|
+ Storer (), /* Write result to 'answer' file. */
|
|
|
+ Verify (); /* Verify the radiosity solution from solver. */
|
|
|
+void SetUp1 (), /* Set up matrix to solve. */
|
|
|
+ SetUp2 (), /* Set up matrix to solve. */
|
|
|
+ Solver (); /* Solve the radiosity matrix. */
|
|
|
+
|
|
|
+main ()
|
|
|
+{
|
|
|
+ int ok; /* Return code temporary storage. */
|
|
|
+
|
|
|
+ /* Get desired number of seconds: */
|
|
|
+ printf ("Enter the number of seconds that is the goal: ");
|
|
|
+ scanf ("%lg", &goal);
|
|
|
+
|
|
|
+ /* Get lower and upper bounds for n from the standard input device: */
|
|
|
+ do {
|
|
|
+ printf ("Enter a lower bound for n: ");
|
|
|
+ scanf ("%d", &n);
|
|
|
+ if (n <= 0)
|
|
|
+ exit(0);
|
|
|
+ ok = Meter (n, &timing, &work);
|
|
|
+ if (timing >= goal)
|
|
|
+ printf ("Must take less than %g seconds. Took %g.\n",
|
|
|
+ goal, timing);
|
|
|
+ } while (!ok || timing >= goal);
|
|
|
+
|
|
|
+ do {
|
|
|
+ printf ("Enter an upper bound for n: ");
|
|
|
+ scanf ("%d", &nupper);
|
|
|
+ if (nupper <= 0)
|
|
|
+ exit(0);
|
|
|
+ ok = Meter (nupper, &timing, &work);
|
|
|
+ if (timing < goal) {
|
|
|
+ printf ("Must take at least %g seconds. Took %g.\n",
|
|
|
+ goal, timing);
|
|
|
+ n = MAX(nupper, n);
|
|
|
+ }
|
|
|
+ } while (!ok || timing < goal);
|
|
|
+
|
|
|
+ /*
|
|
|
+ * While the [n, nupper] interval is larger than 1, bisect it and
|
|
|
+ * pick a half:
|
|
|
+ */
|
|
|
+ while (nupper - n > 1) {
|
|
|
+ mean = (n + nupper) / 2;
|
|
|
+ ok = Meter (mean, &timing, &work);
|
|
|
+ if (timing < goal)
|
|
|
+ n = mean;
|
|
|
+ else
|
|
|
+ nupper = mean;
|
|
|
+ printf ("New interval: [%d,%d]\n", n, nupper);
|
|
|
+ }
|
|
|
+
|
|
|
+ /* Ensure that most recent run was for n, not nupper. */
|
|
|
+ ok = Meter (n, &timing, &work);
|
|
|
+
|
|
|
+ /* Print out final statistics. */
|
|
|
+ What (n, timing, work);
|
|
|
+}
|
|
|
+
|
|
|
+/*****************************************************************************/
|
|
|
+/* This routine should be edited to contain information for your system. */
|
|
|
+/*****************************************************************************/
|
|
|
+What (n, timing, work)
|
|
|
+int n;
|
|
|
+double timing, work;
|
|
|
+{
|
|
|
+ int i;
|
|
|
+ static char *info[] = {
|
|
|
+ "Machine: SUN 4/370GX Processor: SPARC",
|
|
|
+ "Memory: 32 MB # of procs: 1",
|
|
|
+ "Cache: 128 KB # used: 1",
|
|
|
+ "NMAX: 512 Clock: 25 MHz",
|
|
|
+ "Disk: .3GB SCSI+.7GB SMD Node name: amssun2",
|
|
|
+ "OS: SUNOS 4.0.3 Timer: Wall, gettimeofday()",
|
|
|
+ "Language: C Alone: yes",
|
|
|
+ "Compiler: cc Run by: M. Carter",
|
|
|
+ "Options: -O Date: 23 May 1990",
|
|
|
+ NULL
|
|
|
+ };
|
|
|
+
|
|
|
+ printf ("\n");
|
|
|
+ for (i = 0 ; info[i] ; i++)
|
|
|
+ puts (info[i]);
|
|
|
+ printf ("M ops: %-13lg Time: %-.3lf seconds\n",
|
|
|
+ work * 1e-6, timing);
|
|
|
+ printf ("n: %-6d MFLOPS: %-.5lg\n",
|
|
|
+ n, (work / timing) * 1e-6);
|
|
|
+ printf ("Approximate data memory use: %d bytes.\n",
|
|
|
+ 8 * n * n + 120 * n + 800);
|
|
|
+}
|
|
|
+
|
|
|
+/*****************************************************************************/
|
|
|
+/* This routine measures time required on a revised LINPACK-type benchmark, */
|
|
|
+/* including input, matrix generation, solution, and output. */
|
|
|
+/* */
|
|
|
+/* John Gustafson, Diane Rover, Michael Carter, and Stephen Elbert */
|
|
|
+/* Ames Laboratory, 3/18/90 */
|
|
|
+/* */
|
|
|
+/* Calls: Reader Reads the problem description from secondary storage. */
|
|
|
+/* Region Partitions box surface into rectangular regions (patches).*/
|
|
|
+/* SetUp1 Sets up equations from patch geometries-parallel faces. */
|
|
|
+/* SetUp2 Sets up equations from patch geometries-orthogonal faces. */
|
|
|
+/* SetUp3 Sets up equations-row normalization and radiant props. */
|
|
|
+/* Solver Solves the equations by LDL factorization. */
|
|
|
+/* Storer Stores solution (patch radiosities) on secondary storage. */
|
|
|
+/* When Returns wall-clock time, in seconds. */
|
|
|
+/*****************************************************************************/
|
|
|
+
|
|
|
+Meter (npatch, timing, work)
|
|
|
+int npatch; /* In, problem size, here the number of equations. */
|
|
|
+double *timing, /* Out, number of seconds measured. */
|
|
|
+ *work; /* Out, work done, here the number of FLOPs. */
|
|
|
+{
|
|
|
+ static
|
|
|
+ double area[NMAX], /* Areas of patches * 8 * pi. */
|
|
|
+ box[7], /* Dimensions of box in x, y, z directions. */
|
|
|
+ coeff[NMAX][NMAX], /* The coefficients of the eqns to solve. */
|
|
|
+ diag[3][NMAX], /* Diag terms of the eqns to solve. (RGB) */
|
|
|
+ emiss[6][3], /* (RGB) emissivities of patches. */
|
|
|
+ place[3][NMAX], /* Width-height-depth position of patches. */
|
|
|
+ result[3][NMAX], /* Answer radiosities (RGB). */
|
|
|
+ rho[6][3], /* (RGB) Reflectivities of patches. */
|
|
|
+ rhs[3][NMAX], /* Right-hand sides of eqns to solve (RGB). */
|
|
|
+ size[2][NMAX]; /* Width-height sizes of patches. */
|
|
|
+ double ops[8], /* Floating-point operation counts. */
|
|
|
+ p[6], /* Number of patches in faces. */
|
|
|
+ sec[8], /* Times for routines, in seconds. */
|
|
|
+ tmp1, tmp2; /* Double temporary variables. */
|
|
|
+ int i, /* Loop counter. */
|
|
|
+ itmp1, /* Integer temporary variable. */
|
|
|
+ non0; /* Index of first nonzero off-diagonal elem. */
|
|
|
+ static
|
|
|
+ int loop[6][2]; /* Patch number ranges for faces. */
|
|
|
+ static char *tasks[] = { /* Names of all the functions in benchmark. */
|
|
|
+ "Reader", "Region",
|
|
|
+ "SetUp1", "SetUp2",
|
|
|
+ "SetUp3", "Solver",
|
|
|
+ "Storer"
|
|
|
+ };
|
|
|
+ static char *format = /* Output line format. */
|
|
|
+ "%6.6s%8.3f%17.0f%14.6f%10.1f %%\n";
|
|
|
+
|
|
|
+ /* First check that npatch lies between 6 and NMAX: */
|
|
|
+ if (npatch < 6) {
|
|
|
+ printf ("Must be at least 6, the number of faces.\n");
|
|
|
+ return (FALSE);
|
|
|
+ }
|
|
|
+ else if (npatch > NMAX) {
|
|
|
+ printf ("Exceeds %d = maximum for this system.\n", NMAX);
|
|
|
+ return (FALSE);
|
|
|
+ }
|
|
|
+
|
|
|
+ /* Ensure that previous 'answer' file is deleted: */
|
|
|
+ unlink ("answer");
|
|
|
+
|
|
|
+ /* Time the tasks, individually and collectively. */
|
|
|
+ sec[0] = When();
|
|
|
+ if (!Reader (box, rho, emiss))
|
|
|
+ return (FALSE);
|
|
|
+ sec[1] = When();
|
|
|
+ if (!Region (npatch, loop, box, place, size, area))
|
|
|
+ return (FALSE);
|
|
|
+ sec[2] = When();
|
|
|
+ SetUp1 (npatch, loop, coeff, place, size);
|
|
|
+ sec[3] = When();
|
|
|
+ SetUp2 (npatch, loop, coeff, place, size);
|
|
|
+ sec[4] = When();
|
|
|
+ if (!SetUp3 (npatch, loop, area, rho, emiss, coeff, diag, rhs))
|
|
|
+ return (FALSE);
|
|
|
+ sec[5] = When();
|
|
|
+ non0 = loop[1][0];
|
|
|
+ Solver (npatch, non0, coeff, diag, rhs, result);
|
|
|
+ sec[6] = When();
|
|
|
+ Storer (npatch, loop, place, size, result);
|
|
|
+ sec[7] = When();
|
|
|
+ *timing = sec[7] - sec[0];
|
|
|
+ for (i = 0 ; i < 7 ; i++)
|
|
|
+ sec[i] = sec[i+1] - sec[i];
|
|
|
+
|
|
|
+ /* Assess floating-point work done by each routine called, and total: */
|
|
|
+ /* Note the ops counts are talleyed into a double array, and there */
|
|
|
+ /* some strange casts to double in some equations. This is to */
|
|
|
+ /* prevent integer overflow. */
|
|
|
+ itmp1 = 0;
|
|
|
+ tmp1 = 0.0;
|
|
|
+ for (i = 0 ; i < 6 ; i++) {
|
|
|
+ p[i] = loop[i][1] - loop[i][0] + 1;
|
|
|
+ tmp1 += p[i] * p[i];
|
|
|
+ itmp1 += sqrt(p[i] * box[i] / box[i + 1]) + 0.5;
|
|
|
+ }
|
|
|
+ tmp2 = p[0] * p[3] + p[1] * p[4] + p[2] * p[5];
|
|
|
+ ops[0] = 258;
|
|
|
+ ops[1] = 154 + (double) 8 * itmp1 + npatch;
|
|
|
+ ops[2] = 6 + 532 * tmp2;
|
|
|
+ ops[3] = 8*npatch + 370 * ((double) npatch * npatch - tmp1 - 2*tmp2) / 2.0;
|
|
|
+ ops[4] = 72 + (double) 9 * npatch + (double) npatch * npatch - tmp1;
|
|
|
+ ops[5] = npatch * (npatch * ((double) npatch + 7.5) - 2.5) - 21
|
|
|
+ + (non0+1) * ((non0+1) * (2 * ((double) non0+1) - 16.5) + 35.5)
|
|
|
+ + (non0+1) * npatch * (9 - 3 * ((double) non0+1));
|
|
|
+ ops[6] = 48 * npatch;
|
|
|
+ *work = ops[0] + ops[1] + ops[2] + ops[3] + ops[4] + ops[5] + ops[6];
|
|
|
+
|
|
|
+ /* Display timing-work-speed breakdown by routine. */
|
|
|
+ printf ("%d patches:\n", npatch);
|
|
|
+ printf (" Task Seconds Operations MFLOPS %% of Time\n");
|
|
|
+ for (i = 0 ; i < 7 ; i++) {
|
|
|
+ if (sec[i] == 0.0)
|
|
|
+ sec[i] = 0.001;
|
|
|
+ printf (format, tasks[i], sec[i], ops[i], (ops[i] / sec[i]) * 1e-6,
|
|
|
+ 100.0 * sec[i] / *timing);
|
|
|
+ }
|
|
|
+ printf (format, "TOTALS", *timing, *work, (*work / *timing) * 1e-6, 100.0);
|
|
|
+ Verify (npatch, coeff, diag, rhs, result);
|
|
|
+
|
|
|
+ return (TRUE);
|
|
|
+}
|
|
|
+
|
|
|
+/*****************************************************************************/
|
|
|
+/* This function should return the actual, wall clock time (not CPU time) */
|
|
|
+/* in seconds as accurately as possible. Change it to your system timer. */
|
|
|
+/*****************************************************************************/
|
|
|
+double
|
|
|
+When()
|
|
|
+{
|
|
|
+ struct timeval tp;
|
|
|
+ struct timezone tzp;
|
|
|
+ gettimeofday (&tp, &tzp);
|
|
|
+ return ((double) tp.tv_sec + (double) tp.tv_usec * 1e-6);
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+/*****************************************************************************/
|
|
|
+/* The following routine reads in the problem description from secondary */
|
|
|
+/* storage, and checks that numbers are in reasonable ranges. */
|
|
|
+/*****************************************************************************/
|
|
|
+Reader (box, rho, emiss)
|
|
|
+double box[], /* Out: Dimensions of box in x, y, z directions. */
|
|
|
+ rho[][3], /* Out: (RGB) Reflectivities of patches. */
|
|
|
+ emiss[][3]; /* Out: (RGB) emissivities of patches. */
|
|
|
+{
|
|
|
+ /*
|
|
|
+ * Local variables:
|
|
|
+ * infile Device number for input file.
|
|
|
+ * i, j Loop counters.
|
|
|
+ * tmp1 Maximum emissivity, to check that emissivities are not all 0.
|
|
|
+ */
|
|
|
+ int i, j, /* Loop counters. */
|
|
|
+ n; /* Number of args fscanf()'ed from file. */
|
|
|
+ double tmp1; /* Maximum emissivity. */
|
|
|
+ FILE *infile; /* Input file pointer. */
|
|
|
+ char buff[81]; /* Buffer used to eat a line of input. */
|
|
|
+
|
|
|
+ /* Open the input file and read in the data. */
|
|
|
+ if ((infile = fopen ("geom", "r")) == NULL) {
|
|
|
+ printf ("slalom: 'geom' geometry file not found.\n");
|
|
|
+ exit (1);
|
|
|
+ }
|
|
|
+
|
|
|
+ /* Read the box coordinates and error check. */
|
|
|
+ n = 0;
|
|
|
+ for (i = 0 ; i < 3 ; i++) {
|
|
|
+ n += fscanf (infile, "%lg", &box[i]);
|
|
|
+ }
|
|
|
+ fgets (buff, 80, infile); /* Eat the rest of the line. */
|
|
|
+ if (n != 3) {
|
|
|
+ printf ("Must specify exactly 3 box coordinates.\n");
|
|
|
+ exit(1);
|
|
|
+ }
|
|
|
+
|
|
|
+ /* Read the reflectivities and error check. */
|
|
|
+ n = 0;
|
|
|
+ for (j = 0 ; j < 3 ; j++) {
|
|
|
+ for (i = 0 ; i < 6 ; i++) {
|
|
|
+ n += fscanf (infile, "%lg", &rho[i][j]);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ fgets (buff, 80, infile); /* Eat the rest of the line. */
|
|
|
+ if (n != 18) {
|
|
|
+ printf ("Must specify exactly 18 box coordinates.\n");
|
|
|
+ exit(1);
|
|
|
+ }
|
|
|
+
|
|
|
+ /* Read the emissivities and error check. */
|
|
|
+ n = 0;
|
|
|
+ for (j = 0 ; j < 3 ; j++) {
|
|
|
+ for (i = 0 ; i < 6 ; i++) {
|
|
|
+ n += fscanf (infile, "%lg", &emiss[i][j]);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ fgets (buff, 80, infile); /* Eat the rest of the line. */
|
|
|
+ if (n != 18) {
|
|
|
+ printf ("Must specify exactly 18 box coordinates.\n");
|
|
|
+ exit(1);
|
|
|
+ }
|
|
|
+ fclose (infile);
|
|
|
+
|
|
|
+ /* Now sanity-check the values that were just read. */
|
|
|
+ for (j = 0 ; j < 3 ; j++) {
|
|
|
+ if (box[j] < 1.0 || box[j] >= 100.0) {
|
|
|
+ printf ("Box dimensions must be between 1 and 100.\n");
|
|
|
+ return (FALSE);
|
|
|
+ }
|
|
|
+ box[j+3] = box[j];
|
|
|
+
|
|
|
+ tmp1 = 0.0;
|
|
|
+ for (i = 0 ; i < 6 ; i++) {
|
|
|
+ if (rho[i][j] < 0.000 || rho[i][j] > 0.999) {
|
|
|
+ printf ("Reflectivities must be between .000 and .999.\n");
|
|
|
+ return (FALSE);
|
|
|
+ }
|
|
|
+ if (emiss[i][j] < 0.0) {
|
|
|
+ printf ("Emissivity cannot be negative.\n");
|
|
|
+ return (FALSE);
|
|
|
+ }
|
|
|
+ if (tmp1 < emiss[i][j])
|
|
|
+ tmp1 = emiss[i][j];
|
|
|
+ }
|
|
|
+ if (tmp1 == 0.0) {
|
|
|
+ printf ("Emissivities are zero. Problem is trivial.\n");
|
|
|
+ return (FALSE);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ box[6] = box[3];
|
|
|
+ return (TRUE);
|
|
|
+}
|
|
|
+
|
|
|
+/*****************************************************************************/
|
|
|
+/* The following routine decomposes the surface of a variable-sized box */
|
|
|
+/* into patches that are as nearly equal in size and square as possible. */
|
|
|
+/*****************************************************************************/
|
|
|
+Region (npatch, loop, box, place, size, area)
|
|
|
+int npatch, /* In: Problem size. */
|
|
|
+ loop[][2]; /* Out: Patch number ranges for faces. */
|
|
|
+double area[], /* Out: 8pi * areas of the patches. */
|
|
|
+ box[], /* In: Dimensions of box in x, y, z directions. */
|
|
|
+ place[][NMAX], /* Out: Width-height-depth positions of patches. */
|
|
|
+ size[][NMAX]; /* Out: Width-height sizes of patches. */
|
|
|
+{
|
|
|
+
|
|
|
+
|
|
|
+ int icol, /* Loop counter over the number of columns. */
|
|
|
+ ipatch, /* Loop counter over the number of patches. */
|
|
|
+ iface, /* Loop counter over the number of faces. */
|
|
|
+ itmp1, /* Integer temporary variables. */
|
|
|
+ itmp2, /* Integer temporary variables. */
|
|
|
+ last, /* Inner loop ending value. */
|
|
|
+ lead, /* Inner loop starting value. */
|
|
|
+ numcol, /* Number of columns on faces. */
|
|
|
+ numpat, /* Number of patches on a face. */
|
|
|
+ numrow; /* Number of rows of patches in a column. */
|
|
|
+ double height, /* Height of a patch within a column. */
|
|
|
+ tmp1, /* double temporary variables. */
|
|
|
+ tmp2, /* double temporary variables. */
|
|
|
+ tmp3, /* double temporary variables. */
|
|
|
+ tmp4, /* double temporary variables. */
|
|
|
+ width; /* Width of a column of patches. */
|
|
|
+
|
|
|
+ /* Allocate patches to each face, proportionate to area of each face. */
|
|
|
+ tmp1 = 2.0 * (box[0] * box[1] + box[1] * box[2] + box[2] * box[0]);
|
|
|
+ tmp2 = 0.0;
|
|
|
+ tmp3 = npatch;
|
|
|
+ loop[0][0] = 0;
|
|
|
+ for (iface = 0 ; iface < 5 ; iface++) {
|
|
|
+ tmp2 = tmp2 + box[iface] * box[iface + 1];
|
|
|
+ loop[iface][1] = (int) (tmp3 * tmp2 / tmp1 + 0.5) - 1;
|
|
|
+ loop[iface + 1][0] = loop[iface][1] + 1;
|
|
|
+ }
|
|
|
+ loop[5][1] = npatch - 1;
|
|
|
+
|
|
|
+ /* Subdivide each face into numpat patches. */
|
|
|
+ for (iface = 0 ; iface < 6 ; iface++) {
|
|
|
+ numpat = loop[iface][1] - loop[iface][0] + 1;
|
|
|
+ tmp3 = 0.0;
|
|
|
+ if (iface >= 3)
|
|
|
+ tmp3 = box[iface-1];
|
|
|
+ numcol = (int) (sqrt(numpat * box[iface] / box[iface + 1]) + 0.5);
|
|
|
+ if (numcol > numpat)
|
|
|
+ numcol = numpat;
|
|
|
+ if (numcol == 0)
|
|
|
+ numcol = 1;
|
|
|
+ width = box[iface] / numcol;
|
|
|
+ itmp1 = numcol - 1;
|
|
|
+ tmp1 = 0.0;
|
|
|
+ for (icol = 0 ; icol < numcol ; icol++) {
|
|
|
+ itmp2 = itmp1 / numcol;
|
|
|
+ numrow = (itmp1 + numpat) / numcol - itmp2;
|
|
|
+ if (numrow == 0) {
|
|
|
+ printf ("Eccentric box requires more patches.\n");
|
|
|
+ return (FALSE);
|
|
|
+ }
|
|
|
+ height = box[iface + 1] / numrow;
|
|
|
+ tmp2 = 0.0;
|
|
|
+ tmp4 = width * height * (8.0 * M_PI);
|
|
|
+ lead = loop[iface][0] + itmp2;
|
|
|
+ last = lead + numrow;
|
|
|
+ for (ipatch = lead ; ipatch < last ; ipatch++) {
|
|
|
+ size[0][ipatch] = width;
|
|
|
+ size[1][ipatch] = height;
|
|
|
+ place[0][ipatch] = tmp1;
|
|
|
+ place[1][ipatch] = tmp2;
|
|
|
+ place[2][ipatch] = tmp3;
|
|
|
+ area[ipatch] = tmp4;
|
|
|
+ tmp2 = tmp2 + height;
|
|
|
+ }
|
|
|
+ tmp1 = tmp1 + width;
|
|
|
+ itmp1 = itmp1 + numpat;
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ return (TRUE);
|
|
|
+}
|
|
|
+
|
|
|
+/*****************************************************************************/
|
|
|
+/* This routine sets up the radiosity matrix for parallel patches. */
|
|
|
+/*****************************************************************************/
|
|
|
+void
|
|
|
+SetUp1 (npatch, loop, coeff, place, size)
|
|
|
+int npatch, /* In: Problem size. */
|
|
|
+ loop[][2]; /* In: Patch number ranges for faces. */
|
|
|
+double coeff[][NMAX], /* Out: The coefficients of the eqns to solve. */
|
|
|
+ place[][NMAX], /* In: Width-height-depth positions of patches. */
|
|
|
+ size[][NMAX]; /* In: Width-height sizes of patches. */
|
|
|
+{
|
|
|
+ int i, j, k, /* General loop counters. */
|
|
|
+ m, n, /* General loop counters. */
|
|
|
+ iface, /* Loop counter over the number of faces. */
|
|
|
+ ipatch, /* Loop counter over the number of patches. */
|
|
|
+ jface, /* Face coupled to iface when computing mat. elems. */
|
|
|
+ jpatch; /* Patch coupled to ipatch when computing mat. elems.*/
|
|
|
+ double d[2][2][2], /* Point-to-point couplings between patch corners. */
|
|
|
+ d2[2][2][2],/* Squares of d values, to save recomputation. */
|
|
|
+ tmp1, tmp2, /* Double temporary variables. */
|
|
|
+ tmp3, tmp4, /* Double temporary variables. */
|
|
|
+ tmp5, tmp6, /* Double temporary variables. */
|
|
|
+ tmp7, tmp8; /* Double temporary variables. */
|
|
|
+
|
|
|
+ for (iface = 0 ; iface < 3 ; iface++) {
|
|
|
+ jface = iface + 3;
|
|
|
+ tmp1 = place[2][loop[jface][0]] * place[2][loop[jface][0]];
|
|
|
+ tmp6 = tmp1 + tmp1;
|
|
|
+ for (ipatch = loop[iface][0] ; ipatch <= loop[iface][1] ; ipatch++) {
|
|
|
+ for (jpatch=loop[jface][0] ; jpatch <= loop[jface][1] ; jpatch++) {
|
|
|
+ for (j = 0 ; j < 2 ; j++) {
|
|
|
+ d [0][0][j] = place[j][jpatch] - place[j][ipatch];
|
|
|
+ d [1][0][j] = d[0][0][j] + size[j][jpatch];
|
|
|
+ d [0][1][j] = d[0][0][j] - size[j][ipatch];
|
|
|
+ d [1][1][j] = d[1][0][j] - size[j][ipatch];
|
|
|
+ d2[0][0][j] = d[0][0][j] * d[0][0][j];
|
|
|
+ d2[1][0][j] = d[1][0][j] * d[1][0][j];
|
|
|
+ d2[0][1][j] = d[0][1][j] * d[0][1][j];
|
|
|
+ d2[1][1][j] = d[1][1][j] * d[1][1][j];
|
|
|
+ }
|
|
|
+ tmp2 = 0.0;
|
|
|
+ for (m = 0 ; m < 2 ; m++) {
|
|
|
+ for (i = 0 ; i < 2 ; i++) {
|
|
|
+ tmp3 = d2[m][i][1] + tmp1;
|
|
|
+ tmp4 = sqrt(tmp3);
|
|
|
+ tmp5 = 1.0 / tmp4;
|
|
|
+ tmp8 = 0.0;
|
|
|
+ for (k = 0 ; k < 2 ; k++) {
|
|
|
+ for (n = 0 ; n < 2 ; n++) {
|
|
|
+ tmp7 = d[k][n][0];
|
|
|
+ tmp8 = -tmp7 * atan(tmp7 * tmp5) - tmp8;
|
|
|
+ }
|
|
|
+ tmp8 = -tmp8;
|
|
|
+ }
|
|
|
+ tmp2 = -4.0 * tmp4 * tmp8 - tmp2 - tmp6 *
|
|
|
+ log(((d2[1][0][0] + tmp3) * (d2[0][1][0] + tmp3)) /
|
|
|
+ ((d2[0][0][0] + tmp3) * (d2[1][1][0] + tmp3)));
|
|
|
+ }
|
|
|
+ tmp2 = -tmp2;
|
|
|
+ }
|
|
|
+ for (m = 0 ; m < 2 ; m++) {
|
|
|
+ for (i = 0 ; i < 2 ; i++) {
|
|
|
+ tmp4 = sqrt(d2[m][i][0] + tmp1);
|
|
|
+ tmp5 = 1.0 / tmp4;
|
|
|
+ tmp8 = 0.0;
|
|
|
+ for (k = 0 ; k < 2 ; k++) {
|
|
|
+ for (n = 0 ; n < 2 ; n++) {
|
|
|
+ tmp7 = d[k][n][1];
|
|
|
+ tmp8 = -tmp7 * atan(tmp7 * tmp5) - tmp8;
|
|
|
+ }
|
|
|
+ tmp8 = -tmp8;
|
|
|
+ }
|
|
|
+ tmp2 = -4.0 * tmp4 * tmp8 - tmp2;
|
|
|
+ }
|
|
|
+ tmp2 = -tmp2;
|
|
|
+ }
|
|
|
+ coeff[ipatch][jpatch] = tmp2;
|
|
|
+ coeff[jpatch][ipatch] = tmp2;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+/*****************************************************************************/
|
|
|
+/* This routine sets up the radiosity matrix for orthogonal patches. */
|
|
|
+/*****************************************************************************/
|
|
|
+void
|
|
|
+SetUp2 (npatch, loop, coeff, place, size)
|
|
|
+int npatch, /* In: Problem size. */
|
|
|
+ loop[][2]; /* In: Patch number ranges for faces. */
|
|
|
+double coeff[][NMAX], /* Out: The coefficients of the eqns to solve. */
|
|
|
+ place[][NMAX], /* In: Width-height-depth positions of patches. */
|
|
|
+ size[][NMAX]; /* In: Width-height sizes of patches. */
|
|
|
+{
|
|
|
+ int m, /* General loop counters. */
|
|
|
+ iface, /* Loop counter over the number of faces. */
|
|
|
+ ipatch, /* Loop counter over the number of patches. */
|
|
|
+ jface, /* Face coupled to iface when computing mat. elems. */
|
|
|
+ jpatch; /* Patch coupled to ipatch when computing mat. elems.*/
|
|
|
+
|
|
|
+ double tmpb, tmpa,
|
|
|
+ c11d, c12d, c21d, c22d, c11s, c12s, c21s, c22s,
|
|
|
+ d11d, d12d, d21d, d22d, d11s, d12s, d21s, d22s,
|
|
|
+ d11i, d12i, d21i, d22i, a10s, a20s, b01s, b02s,
|
|
|
+ e1111, e1211, e2111, e2211, e1112, e1212, e2112, e2212,
|
|
|
+ e1121, e1221, e2121, e2221, e1122, e1222, e2122, e2222;
|
|
|
+
|
|
|
+ for (iface = 0 ; iface < 6 ; iface++) {
|
|
|
+ for (m = 0 ; m < 2 ; m++) {
|
|
|
+ jface = (iface + m + 1) % 6;
|
|
|
+ for (ipatch=loop[iface][0] ; ipatch <= loop[iface][1] ; ipatch++) {
|
|
|
+ a10s = place[m][ipatch] - place[2][loop[jface][0]];
|
|
|
+ a20s = a10s + size[m][ipatch];
|
|
|
+ a10s = a10s * a10s;
|
|
|
+ a20s = a20s * a20s;
|
|
|
+ for (jpatch=loop[jface][0] ; jpatch<=loop[jface][1];jpatch++) {
|
|
|
+ c11d = place[m][jpatch] - place[1-m][ipatch];
|
|
|
+ c12d = c11d + size[m][jpatch];
|
|
|
+ c21d = c11d - size[1-m][ipatch];
|
|
|
+ c22d = c12d - size[1-m][ipatch];
|
|
|
+ c11s = c11d * c11d;
|
|
|
+ c12s = c12d * c12d;
|
|
|
+ c21s = c21d * c21d;
|
|
|
+ c22s = c22d * c22d;
|
|
|
+ b01s = place[1 - m][jpatch] - place[2][ipatch];
|
|
|
+ b02s = b01s + size[1 - m][jpatch];
|
|
|
+
|
|
|
+ /**/
|
|
|
+ /* Bump the term by a small real to avoid
|
|
|
+ /* singularities in coupling function:
|
|
|
+ /**/
|
|
|
+ b01s = b01s * b01s + 1e-35;
|
|
|
+ b02s = b02s * b02s + 1e-35;
|
|
|
+ d11s = a10s + b01s;
|
|
|
+ d12s = a10s + b02s;
|
|
|
+ d21s = a20s + b01s;
|
|
|
+ d22s = a20s + b02s;
|
|
|
+ d11d = sqrt(d11s);
|
|
|
+ d12d = sqrt(d12s);
|
|
|
+ d21d = sqrt(d21s);
|
|
|
+ d22d = sqrt(d22s);
|
|
|
+ d11i = 1.0 / d11d;
|
|
|
+ d12i = 1.0 / d12d;
|
|
|
+ d21i = 1.0 / d21d;
|
|
|
+ d22i = 1.0 / d22d;
|
|
|
+
|
|
|
+ tmpa = d11d * ( c11d * atan (c11d * d11i)
|
|
|
+ - c12d * atan (c12d * d11i)
|
|
|
+ - c21d * atan (c21d * d11i)
|
|
|
+ + c22d * atan (c22d * d11i))
|
|
|
+ + d12d * (-c11d * atan (c11d * d12i)
|
|
|
+ + c12d * atan (c12d * d12i)
|
|
|
+ + c21d * atan (c21d * d12i)
|
|
|
+ - c22d * atan (c22d * d12i))
|
|
|
+ + d21d * (-c11d * atan (c11d * d21i)
|
|
|
+ + c12d * atan (c12d * d21i)
|
|
|
+ + c21d * atan (c21d * d21i)
|
|
|
+ - c22d * atan (c22d * d21i))
|
|
|
+ + d22d * ( c11d * atan (c11d * d22i)
|
|
|
+ - c12d * atan (c12d * d22i)
|
|
|
+ - c21d * atan (c21d * d22i)
|
|
|
+ + c22d * atan (c22d * d22i));
|
|
|
+
|
|
|
+ e1111 = c11s + d11s;
|
|
|
+ e1211 = c12s + d11s;
|
|
|
+ e2111 = c21s + d11s;
|
|
|
+ e2211 = c22s + d11s;
|
|
|
+ e1112 = c11s + d12s;
|
|
|
+ e1212 = c12s + d12s;
|
|
|
+ e2112 = c21s + d12s;
|
|
|
+ e2212 = c22s + d12s;
|
|
|
+ e1121 = c11s + d21s;
|
|
|
+ e1221 = c12s + d21s;
|
|
|
+ e2121 = c21s + d21s;
|
|
|
+ e2221 = c22s + d21s;
|
|
|
+ e1122 = c11s + d22s;
|
|
|
+ e1222 = c12s + d22s;
|
|
|
+ e2122 = c21s + d22s;
|
|
|
+ e2222 = c22s + d22s;
|
|
|
+
|
|
|
+ tmpb = c11s * log( e1111 * e1122 / (e1112 * e1121))
|
|
|
+ - c12s * log( e1211 * e1222 / (e1212 * e1221))
|
|
|
+ - c21s * log( e2111 * e2122 / (e2112 * e2121))
|
|
|
+ + c22s * log( e2211 * e2222 / (e2212 * e2221))
|
|
|
+ - d11s * log( e1111 * e2211 / (e1211 * e2111))
|
|
|
+ + d12s * log( e1112 * e2212 / (e1212 * e2112))
|
|
|
+ + d21s * log( e1121 * e2221 / (e1221 * e2121))
|
|
|
+ - d22s * log( e1122 * e2222 / (e1222 * e2122));
|
|
|
+
|
|
|
+ coeff[ipatch][jpatch] = fabs(4.0 * tmpa + tmpb);
|
|
|
+ coeff[jpatch][ipatch] = coeff[ipatch][jpatch];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+/*****************************************************************************/
|
|
|
+/* This routine sets up the radiosity matrix... normalizes row sums to 1, */
|
|
|
+/* and includes terms derived from reflectivites and emissivities of faces. */
|
|
|
+/*****************************************************************************/
|
|
|
+SetUp3 (npatch, loop, area, rho, emiss, coeff, diag, rhs)
|
|
|
+int npatch, /* In: Problem size. */
|
|
|
+ loop[][2]; /* In: Patch number ranges for faces. */
|
|
|
+double area[], /* In: 8 * pi * areas of the patches. */
|
|
|
+ rho[][3], /* In: (RGB) Reflectivities of the face interiors. */
|
|
|
+ emiss[][3], /* In: (RGB) Emissivities of the face interiors. */
|
|
|
+ coeff[][NMAX], /* Out: The coefficients of the eqns to solve. */
|
|
|
+ diag[][NMAX], /* Out: (RGB) Diagonal terms of the system. */
|
|
|
+ rhs[][NMAX]; /* Out: (RGB) Right-hand sides of system to solve. */
|
|
|
+{
|
|
|
+
|
|
|
+ /*
|
|
|
+ * Local variables:
|
|
|
+ * iface Loop counter over the number of faces.
|
|
|
+ * ipatch Outer loop counter over the number of patches.
|
|
|
+ * j Loop counter over each color (R-G-B).
|
|
|
+ * jpatch Inner loop counter over the number of patches.
|
|
|
+ * tmp1 double temporary variable.
|
|
|
+ * vtmp1-2 double vector temporary variables.
|
|
|
+ */
|
|
|
+ int j, /* (RGB) Loop counter over each color. */
|
|
|
+ iface, /* Loop counter over the number of faces. */
|
|
|
+ ipatch, /* Outer loop counter over the number of patches. */
|
|
|
+ jpatch; /* Inner loop counter over the number of patches. */
|
|
|
+ double tmp1, /* Double temporary variable. */
|
|
|
+ vtmp1[3], /* Double vector temporary variables. */
|
|
|
+ vtmp2[3]; /* Double vector temporary variables. */
|
|
|
+
|
|
|
+ /* Ensure that row sums to 1, and put in reflectivities (rho) and */
|
|
|
+ /* emissivities. */
|
|
|
+ for (iface = 0 ; iface < 6 ; iface++) {
|
|
|
+ for (j = 0 ; j < 3 ; j++) {
|
|
|
+ vtmp1[j] = 1.0 / rho[iface][j];
|
|
|
+ vtmp2[j] = emiss[iface][j] * vtmp1[j];
|
|
|
+ }
|
|
|
+ for (ipatch = loop[iface][0] ; ipatch <= loop[iface][1] ; ipatch++) {
|
|
|
+ tmp1 = 0.0;
|
|
|
+ for (jpatch = 0 ; jpatch < loop[iface][0] ; jpatch++) {
|
|
|
+ tmp1 += coeff[ipatch][jpatch];
|
|
|
+ }
|
|
|
+ for (jpatch = loop[iface][1]+1 ; jpatch < npatch ; jpatch++) {
|
|
|
+ tmp1 += coeff[ipatch][jpatch];
|
|
|
+ }
|
|
|
+ /* Make sure row sum (total form factor) is close to 1: */
|
|
|
+ if (fabs(tmp1 - area[ipatch]) > (0.5e-9 * tmp1)) {
|
|
|
+ printf ("Total form factor is too far from unity.\n");
|
|
|
+ return (FALSE);
|
|
|
+ }
|
|
|
+ tmp1 = -tmp1;
|
|
|
+ /* Set coplanar patch interactions to zero. */
|
|
|
+ for (jpatch=loop[iface][0] ; jpatch <= loop[iface][1] ; jpatch++) {
|
|
|
+ coeff[ipatch][jpatch] = 0.0;
|
|
|
+ }
|
|
|
+ /* Assign diagonal entries and right-hand sides. */
|
|
|
+ for (j = 0 ; j < 3 ; j++) {
|
|
|
+ diag[j][ipatch] = vtmp1[j] * tmp1;
|
|
|
+ rhs[j][ipatch] = vtmp2[j] * tmp1;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return (TRUE);
|
|
|
+}
|
|
|
+
|
|
|
+/*****************************************************************************/
|
|
|
+/* This routine factors and backsolves a real, symmetric, near-dense matrix */
|
|
|
+/* by LDL factorization. No pivoting; the matrix is diagonally dominant. */
|
|
|
+/*****************************************************************************/
|
|
|
+void
|
|
|
+Solver (npatch, non0, coeff, diag, rhs, result)
|
|
|
+int npatch, /* In: Problem size. */
|
|
|
+ non0; /* In: Index of first nonzero off-diagonal mat. elem.*/
|
|
|
+double coeff[][NMAX], /* In/Out: The coefficients of the eqns to solve. */
|
|
|
+ diag[][NMAX], /* Out: (RGB) Diagonal terms of the system. */
|
|
|
+ rhs[][NMAX], /* In: (RGB) Right-hand sides of system to solve. */
|
|
|
+ result[][NMAX]; /* Out: (RGB) solution radiosities. */
|
|
|
+{
|
|
|
+ int i, j, /* General loop counters. */
|
|
|
+ k, m; /* General loop counters. */
|
|
|
+ double tmp1; /* Double temporary variable. */
|
|
|
+
|
|
|
+ /* Load lower triangle of coefficients, diagonal, and solution vector. */
|
|
|
+ for (m = 0 ; m < 3 ; m++) {
|
|
|
+ for (i = non0 ; i < npatch ; i++) {
|
|
|
+ coeff[i][i] = diag[m][i];
|
|
|
+ result[m][i] = rhs[m][i];
|
|
|
+ for (j = 0 ; j < i ; j++) {
|
|
|
+ coeff[i][j] = coeff[j][i];
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ /* Factor matrix, writing factors on top of original matrix. */
|
|
|
+ for (j = 0 ; j < non0 ; j++) {
|
|
|
+ coeff[j][j] = 1.0 / diag[m][j];
|
|
|
+ result[m][j] = rhs[m][j];
|
|
|
+ }
|
|
|
+
|
|
|
+ for (j = non0 ; j < npatch ; j++) {
|
|
|
+ for (k = non0 ; k < j ; k++) {
|
|
|
+ coeff[j][k] -= Ddot (k, &coeff[k][0], 1, &coeff[j][0], 1);
|
|
|
+ }
|
|
|
+ for (k = 0 ; k < j ; k++) {
|
|
|
+ tmp1 = coeff[j][k];
|
|
|
+ coeff[j][k] = tmp1 * coeff[k][k];
|
|
|
+ coeff[j][j] -= tmp1 * coeff[j][k];
|
|
|
+ }
|
|
|
+ coeff[j][j] = 1.0 / coeff[j][j];
|
|
|
+ }
|
|
|
+
|
|
|
+ /* Backsolve, in three stages (for L, D, and L transpose). */
|
|
|
+ for (k = non0 ; k < npatch ; k++) {
|
|
|
+ result[m][k] -= Ddot (k, &result[m][0], 1, &coeff[k][0], 1);
|
|
|
+ }
|
|
|
+
|
|
|
+ for (k = 0 ; k < npatch ; k++) {
|
|
|
+ result[m][k] *= coeff[k][k];
|
|
|
+ }
|
|
|
+
|
|
|
+ for (k = npatch - 2 ; k >= non0 ; k--) {
|
|
|
+ result[m][k] -= Ddot (npatch-(k+1), &result[m][k+1], 1,
|
|
|
+ &coeff[k+1][k], NMAX);
|
|
|
+ }
|
|
|
+
|
|
|
+ for (k = non0 - 1 ; k >= 0 ; k--) {
|
|
|
+ result[m][k] -= Ddot (npatch-non0, &result[m][non0], 1,
|
|
|
+ &coeff[non0][k], NMAX);
|
|
|
+ }
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+/*****************************************************************************/
|
|
|
+/* The following routine writes the answer to secondary storage. */
|
|
|
+/*****************************************************************************/
|
|
|
+Storer (npatch, loop, place, size, result)
|
|
|
+int npatch, /* In: Problem size. */
|
|
|
+ loop[][2]; /* In: Patch number ranges for faces. */
|
|
|
+double result[][NMAX], /* In: (RGB) Radiosity solutions. */
|
|
|
+ place[][NMAX], /* In: Width-height-depth positions of patches. */
|
|
|
+ size[][NMAX]; /* In: Width-height sizes of patches. */
|
|
|
+{
|
|
|
+ int i, /* General loop counter. */
|
|
|
+ iface, /* Loop counter over number of faces. */
|
|
|
+ ipatch; /* Loop counter of number of patches within a face. */
|
|
|
+ FILE *outfile; /* Output file pointer. */
|
|
|
+
|
|
|
+ /* Write patch geometry to 'answer' file. */
|
|
|
+ if ((outfile = fopen("answer", "w")) == NULL) {
|
|
|
+ printf ("Unable to open 'answer' file.\n");
|
|
|
+ exit (1);
|
|
|
+ }
|
|
|
+ fprintf (outfile, "%d patches:\n", npatch);
|
|
|
+ fprintf (outfile,
|
|
|
+ " Patch Face Position in w, h, d Width Height\n");
|
|
|
+ for (iface = 0 ; iface < 6 ; iface++) {
|
|
|
+ for (ipatch = loop[iface][0] ; ipatch <= loop[iface][1] ; ipatch++) {
|
|
|
+ fprintf (outfile,
|
|
|
+ "%5d %4d%11.5lf%11.5lf%11.5lf %11.5lf%11.5lf\n",
|
|
|
+ ipatch+1, iface+1,
|
|
|
+ place[0][ipatch],
|
|
|
+ place[1][ipatch],
|
|
|
+ place[2][ipatch],
|
|
|
+ size[0][ipatch],
|
|
|
+ size[1][ipatch]);
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ /* Write patch radiosities to 'answer' file. */
|
|
|
+ fprintf (outfile, "\n Patch Face Radiosities\n");
|
|
|
+ for (iface = 0 ; iface < 6 ; iface++) {
|
|
|
+ for (ipatch = loop[iface][0] ; ipatch <= loop[iface][1] ; ipatch++) {
|
|
|
+ fprintf (outfile, "%5d %4d%12.8lf%12.8lf%12.8lf\n",
|
|
|
+ ipatch+1, iface+1,
|
|
|
+ result[0][ipatch],
|
|
|
+ result[1][ipatch],
|
|
|
+ result[2][ipatch]);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ fclose(outfile);
|
|
|
+}
|
|
|
+
|
|
|
+/*****************************************************************************/
|
|
|
+/* This routine verifies that the computed radiosities satisfy the equations.*/
|
|
|
+/* */
|
|
|
+/* John Gustafson, Diane Rover, Michael Carter, and Stephen Elbert */
|
|
|
+/* Ames Laboratory, 3/18/90 */
|
|
|
+/*****************************************************************************/
|
|
|
+Verify (npatch, coeff, diag, rhs, result)
|
|
|
+int npatch; /* In: Problem size. */
|
|
|
+double coeff[][NMAX], /* In: The coefficients of the eqns to solve. */
|
|
|
+ diag[][NMAX], /* In: (RGB) Diagonal terms of the system. */
|
|
|
+ rhs[][NMAX], /* In: (RGB) Right-hand sides of system to solve. */
|
|
|
+ result[][NMAX]; /* In: (RGB) Radiosity solutions. */
|
|
|
+{
|
|
|
+ double tmp1, tmp2; /* Double temporary variables. */
|
|
|
+ double anorm, /* Norm accumulation variable. */
|
|
|
+ xnorm; /* Norm accumulation variable. */
|
|
|
+ int i, j, m; /* General loop counters. */
|
|
|
+
|
|
|
+ tmp1 = 0.0;
|
|
|
+ for (m = 0 ; m < 3 ; m++) {
|
|
|
+ /* Copy lower triangle of coefficients to upper triangle, */
|
|
|
+ /* and load diagonal. */
|
|
|
+ for (i = 0 ; i < npatch ; i++) {
|
|
|
+ coeff[i][i] = diag[m][i];
|
|
|
+ for (j = 0 ; j < i ; j++) {
|
|
|
+ coeff[i][j] = coeff[j][i];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ /* Multiply matrix by solution vector, and accum. norm of residual. */
|
|
|
+ anorm = xnorm = 0.0;
|
|
|
+ for (j = 0 ; j < npatch ; j++) {
|
|
|
+ tmp2 = rhs[m][j];
|
|
|
+ for (i = 0 ; i < npatch ; i++) {
|
|
|
+ tmp2 -= (coeff[j][i] * result[m][i]);
|
|
|
+ anorm = MAX(anorm, fabs(coeff[j][i]));
|
|
|
+ }
|
|
|
+ xnorm = MAX(xnorm, fabs(result[m][j]));
|
|
|
+ tmp1 += fabs(tmp2);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ /* printf ("anorm = %g xnorm = %g\n", anorm, xnorm); */
|
|
|
+ tmp1 /= (anorm * xnorm);
|
|
|
+ if (tmp1 > 3 * EPS) {
|
|
|
+ printf ("Residual is too large: %lg\n", tmp1);
|
|
|
+ return (FALSE);
|
|
|
+ }
|
|
|
+ return (TRUE);
|
|
|
+}
|
|
|
+
|
|
|
+#ifdef SUN4
|
|
|
+
|
|
|
+/*****************************************************************************/
|
|
|
+/* Double precision dot product specifically written for Sun 4/370. */
|
|
|
+/* By Michael Carter and John Gustafson, May 30, 1990 */
|
|
|
+/* This code unrolls the dot product four ways since that's how many */
|
|
|
+/* registers are available on the SPARC. Other RISC system will require */
|
|
|
+/* something very similar. Also, unit stride is take advantage of in the */
|
|
|
+/* form of special cases. */
|
|
|
+/*****************************************************************************/
|
|
|
+double
|
|
|
+Ddot (n, a, ia, b, ib)
|
|
|
+register
|
|
|
+int n, /* Number of elements in vectors. */
|
|
|
+ ia, /* Stride of a vector in ELEMENTS. */
|
|
|
+ ib; /* Stride of b vector in ELEMENTS. */
|
|
|
+register
|
|
|
+double *a, /* Pointer to first vector. */
|
|
|
+ *b; /* Pointer to second vector. */
|
|
|
+{
|
|
|
+ register double sum0 = 0.0,
|
|
|
+ sum1 = 0.0,
|
|
|
+ sum2 = 0.0,
|
|
|
+ sum3 = 0.0;
|
|
|
+ register int m = n & 3;
|
|
|
+ int t;
|
|
|
+
|
|
|
+ /* The ragged cleanup part. */
|
|
|
+ while (m--) {
|
|
|
+ sum0 += *a * *b;
|
|
|
+ a += ia;
|
|
|
+ b += ib;
|
|
|
+ }
|
|
|
+
|
|
|
+ /* The fast pipelined part */
|
|
|
+ n >>= 2;
|
|
|
+ if (ib == 1 && ia != 1) {
|
|
|
+ t = ia;
|
|
|
+ ia = ib;
|
|
|
+ ib = t;
|
|
|
+ t = (int) a;
|
|
|
+ b = a;
|
|
|
+ a = (double *) t;
|
|
|
+ }
|
|
|
+
|
|
|
+ /* We can optimize if one or more strides are equal to 1. */
|
|
|
+ if (ia == 1) {
|
|
|
+ /* This runs if both strides are 1. */
|
|
|
+ if (ib == 1) {
|
|
|
+ ia <<= 2;
|
|
|
+ ib <<= 2;
|
|
|
+ while (n--) {
|
|
|
+ sum0 += a[0] * b[0];
|
|
|
+ sum1 += a[1] * b[1];
|
|
|
+ sum2 += a[2] * b[2];
|
|
|
+ sum3 += a[3] * b[3];
|
|
|
+ a += ia;
|
|
|
+ b += ib;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ /* This runs if stride of a only is equal to 1. */
|
|
|
+ else {
|
|
|
+ ia <<= 2;
|
|
|
+ while (n--) {
|
|
|
+ sum0 += a[0] * *b;
|
|
|
+ b += ib;
|
|
|
+ sum1 += a[1] * *b;
|
|
|
+ b += ib;
|
|
|
+ sum2 += a[2] * *b;
|
|
|
+ b += ib;
|
|
|
+ sum3 += a[3] * *b;
|
|
|
+ a += ia;
|
|
|
+ b += ib;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ /* This runs for the more general case. */
|
|
|
+ /* This is about .5 MFLOPS slower on Sun 4/370 */
|
|
|
+ else {
|
|
|
+ while (n--) {
|
|
|
+ sum0 += *a * *b;
|
|
|
+ a += ia;
|
|
|
+ b += ib;
|
|
|
+ sum1 += *a * *b;
|
|
|
+ a += ia;
|
|
|
+ b += ib;
|
|
|
+ sum2 += *a * *b;
|
|
|
+ a += ia;
|
|
|
+ b += ib;
|
|
|
+ sum3 += *a * *b;
|
|
|
+ a += ia;
|
|
|
+ b += ib;
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ return (sum0 + sum1 + sum2 + sum3);
|
|
|
+}
|
|
|
+
|
|
|
+#else
|
|
|
+
|
|
|
+/*****************************************************************************/
|
|
|
+/* Generic double-precision dot product. Unrolling will help pipelined */
|
|
|
+/* computers. Modify accordingly. */
|
|
|
+/*****************************************************************************/
|
|
|
+double
|
|
|
+Ddot (n, a, ia, b, ib)
|
|
|
+register
|
|
|
+int n, /* Number of elements in vectors. */
|
|
|
+ ia, /* Stride of a vector in ELEMENTS. */
|
|
|
+ ib; /* Stride of b vector in ELEMENTS. */
|
|
|
+register
|
|
|
+double *a, /* Pointer to first vector. */
|
|
|
+ *b; /* Pointer to second vector. */
|
|
|
+{
|
|
|
+ register double sum = 0.0;
|
|
|
+
|
|
|
+ while (n--) {
|
|
|
+ sum += *a * *b;
|
|
|
+ a += ia;
|
|
|
+ b += ib;
|
|
|
+ }
|
|
|
+ return (sum);
|
|
|
+}
|
|
|
+
|
|
|
+#endif
|
Jayke Meijer