#!/bin/sh
# This is a shar archive.
# The rest of this file is a shell script which will extract:
#
# 4_10a1.c 4_10a2.c 4_10b.c 4_10b.cmp 4_10b0.c 4_10b1.cmp 4_10b2.cmp 4_10b3.cmp 4_10c.c 4_10c.cmp 4_10c1.cmp 4_10c2.cmp 4_10c3.cmp inp1 inp2 inp3 main.c makefile
#
# To extract the files from this shell archive file simply
# create a directory for this file, move the archive file
# to it and enter the command
#
#	sh filename
# 
# The files will be extracted automatically.
# Note: Do not use csh.
#
# Archive created: Mon Jul 30 23:02:25 EDT 1990
#
echo x - 4_10a1.c
sed 's/^X//' > 4_10a1.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
double **mkarray(int rows, int cols)
{
    double **ret = new double *[rows];
    for (int i = 0; i < cols; i++)
	ret[i] = new double[cols];
    return ret;
}
!EOF!
ls -l 4_10a1.c
echo x - 4_10a2.c
sed 's/^X//' > 4_10a2.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
void delarray(double **array, int rows, int cols)
{
    for (int i = 0; i < cols; i++)
	delete array[i];
    delete array;
}
!EOF!
ls -l 4_10a2.c
echo x - 4_10b.c
sed 's/^X//' > 4_10b.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
#ifndef VER2		/* DELETE */
#include "4_10b0.c"	/* DELETE */
#else			/* DELETE */
#include "4_10c.c"	/* DELETE */
#endif			/* DELETE */
// invert A and place the result in I
int matinv(double **A, double **I, int nelem)
{
    // initialize identity matrix I
    for (int i = 0; i < nelem; i++)
	{
	for (int j = 0; j < nelem; j++)
	    I[i][j] = 0.0;

	I[i][i] = 1.0;
	}

    // loop diag from 0 to nelem-1
    for (int diag = 0; diag < nelem; diag++)
	{
	if (!dopivot(A, I, diag, nelem))
	    return 0;

	// divide row diag by A[diag][diag] in A and I
	//     ignore A[diag][0..diag-1] since they
	//     are already 0
	double div = A[diag][diag];
	if (div != 1.0)
	    {
	    A[diag][diag] = 1.0;
	    for (int j = diag+1; j < nelem; j++)
		A[diag][j] /= div;

	    for (j = 0; j < nelem; j++)
		I[diag][j] /= div;
	    }

	// loop i from 0 to diag-1 and
	// from diag+1 to nelem-1
	for (i = 0; i < nelem; i++)
	    {
	    if (i == diag)
		continue;

	    // subtract A[i][diag] * row diag
	    // from row i in A and I
	    double sub = A[i][diag];
	    if (sub != 0.0)
		{
		A[i][diag] = 0.0;
		for (int j = diag+1; j < nelem; j++)
		    A[i][j] -= sub * A[diag][j];

		for (j = 0; j < nelem; j++)
		    I[i][j] -= sub * I[diag][j];
		}
	    }
	}

    return 1;
}
!EOF!
ls -l 4_10b.c
echo x - 4_10b.cmp
sed 's/^X//' > 4_10b.cmp << '!EOF!'
!EOF!
ls -l 4_10b.cmp
echo x - 4_10b0.c
sed 's/^X//' > 4_10b0.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
// Exercise 4.10 (version 1)
// Invert a two dimensional (n X n) array
// using the Gauss-Jordan elimination method.
// Return 1 if succeeded, 0 if singular.
#include <swap.h>

// Find A[i][d] which is non zero, i >= d,
// and swap with A[d][d].
// Return FALSE if there are no non-zero entries.
static int dopivot(double **A, double **I,
    int diag, int nelem)
{
    if (A[diag][diag] != 0.0)
        return 1;

    for (int i = diag+1; i < nelem; i++)
	if (A[i][diag] != 0.0)
	    {
	    // swap row i with row diag in A and I
	    swap(A[diag], A[i]);
	    swap(I[diag], I[i]);
	    break;
	    }

    // if no non-zero entries (i == nelem),
    // then failure
    return i < nelem;
}
!EOF!
ls -l 4_10b0.c
echo x - 4_10b1.cmp
sed 's/^X//' > 4_10b1.cmp << '!EOF!'


A	(2 x 2)


	2	3
	1	2


A again	(2 x 2)


	1	0
	0	1


inverse	(2 x 2)


	2	-3
	-1	2


should be	(2 x 2)


	2	-3
	-1	2
!EOF!
ls -l 4_10b1.cmp
echo x - 4_10b2.cmp
sed 's/^X//' > 4_10b2.cmp << '!EOF!'


A	(3 x 3)


	1	2	3
	0	1	2
	0	0	1


A again	(3 x 3)


	1	0	0
	0	1	0
	0	0	1


inverse	(3 x 3)


	1	-2	1
	0	1	-2
	0	0	1


should be	(3 x 3)


	1	-2	1
	0	1	-2
	0	0	1
!EOF!
ls -l 4_10b2.cmp
echo x - 4_10b3.cmp
sed 's/^X//' > 4_10b3.cmp << '!EOF!'


A	(3 x 3)


	2	1	-2
	0	1	1
	-1	0	0


A again	(3 x 3)


	1	0	0
	0	1	0
	0	0	1


inverse	(3 x 3)


	0	0	-1
	0.333333	0.666667	0.666667
	-0.333333	0.333333	-0.666667


should be	(3 x 3)


	0	0	-1
	0.333333	0.666667	0.666667
	-0.333333	0.333333	-0.666667
!EOF!
ls -l 4_10b3.cmp
echo x - 4_10c.c
sed 's/^X//' > 4_10c.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
// Exercise 4.10 (version 2)
// Invert a two dimensional (n X n) array
// using the Gauss-Jordan elimination method
// with partial pivoting.
#include <swap.h>
#include <abs.h>

// Find the maximum |A[i][diag]|, i >= diag.
// Return FALSE if there are no non-zero entries.
static int dopivot(double **A, double **I,
    int diag, int nelem)
{
    double max = abs(A[diag][diag]);
    int swaprow = diag;

    for (int i = diag+1; i < nelem; i++)
	{
	double a_i = abs(A[i][diag]);
	if (a_i > max)
	    {
	    max = a_i;
	    swaprow = i;
	    }
	}

    // if all are zero, then failure
    if (max == 0.0)
	return 0;

    // swap row swaprow with row diag in A and I
    if (swaprow != diag)
	{
	for (int j = diag; j < nelem; j++)
	    swap(A[diag][j], A[swaprow][j]);
	for (j = 0; j < nelem; j++)
	    swap(I[diag][j], I[swaprow][j]);
	}

    return 1;
}
!EOF!
ls -l 4_10c.c
echo x - 4_10c.cmp
sed 's/^X//' > 4_10c.cmp << '!EOF!'
!EOF!
ls -l 4_10c.cmp
echo x - 4_10c1.cmp
sed 's/^X//' > 4_10c1.cmp << '!EOF!'


A	(2 x 2)


	2	3
	1	2


A again	(2 x 2)


	1	0
	0	1


inverse	(2 x 2)


	2	-3
	-1	2


should be	(2 x 2)


	2	-3
	-1	2
!EOF!
ls -l 4_10c1.cmp
echo x - 4_10c2.cmp
sed 's/^X//' > 4_10c2.cmp << '!EOF!'


A	(3 x 3)


	1	2	3
	0	1	2
	0	0	1


A again	(3 x 3)


	1	0	0
	0	1	0
	0	0	1


inverse	(3 x 3)


	1	-2	1
	0	1	-2
	0	0	1


should be	(3 x 3)


	1	-2	1
	0	1	-2
	0	0	1
!EOF!
ls -l 4_10c2.cmp
echo x - 4_10c3.cmp
sed 's/^X//' > 4_10c3.cmp << '!EOF!'


A	(3 x 3)


	2	1	-2
	0	1	1
	-1	0	0


A again	(3 x 3)


	1	0	0
	0	1	0
	0	0	1


inverse	(3 x 3)


	0	0	-1
	0.333333	0.666667	0.666667
	-0.333333	0.333333	-0.666667


should be	(3 x 3)


	0	0	-1
	0.333333	0.666667	0.666667
	-0.333333	0.333333	-0.666667
!EOF!
ls -l 4_10c3.cmp
echo x - inp1
sed 's/^X//' > inp1 << '!EOF!'
2
2 3
1 2

2 -3
-1 2
!EOF!
ls -l inp1
echo x - inp2
sed 's/^X//' > inp2 << '!EOF!'
3
1 2 3
0 1 2
0 0 1


1 -2 1
0 1 -2
0 0 1
!EOF!
ls -l inp2
echo x - inp3
sed 's/^X//' > inp3 << '!EOF!'
3
2 1 -2
0 1 1
-1 0 0


0 0 -1
0.3333333 0.6666667 0.6666667
-0.3333333 0.3333333 -0.6666667
!EOF!
ls -l inp3
echo x - main.c
sed 's/^X//' > main.c << '!EOF!'
/* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
/* The C++ Answer Book */
/* Tony Hansen */
/* All rights reserved. */
#include <stream.h>

double **mkarray(int rows, int cols);
int matinv(double **A, double **I, int n);

void prarray(double **a, int rows, int cols, char *name)
{
    cout << "\n\n" << name << "\t(" << rows << " x " << cols << ")\n\n\n";
    for (int i = 0; i < rows; i++)
	{
	for (int j = 0; j < cols; j++)
	    {
	    cout << "\t" << a[i][j];
	    }
	cout << "\n";
	}
}

void rdarray(double **a, int rows, int cols)
{
    for (int i = 0; i < rows; i++)
        for (int j = 0; j < cols; j++)
	    cin >> a[i][j];
}

main()
{
    int size;
    cin >> size;
    double **A = mkarray(size, size);
    rdarray(A, size, size);
    prarray(A, size, size, "A");
    double **I = mkarray(size, size);
    double **S = mkarray(size, size);
    rdarray(S, size, size);
    int ret = matinv(A, I, size);
    if (ret == 0)
        cout << "singular matrix!\n";
    else
        {
	prarray(A, size, size, "A again");
	prarray(I, size, size, "inverse");
	prarray(S, size, size, "should be");
	}
    return 0;
}
!EOF!
ls -l main.c
echo x - makefile
sed 's/^X//' > makefile << '!EOF!'
CFLAGS = -g +i
CC= CC -I. -I../../CC

all: 4_10b 4_10c

4_10b: main.o 4_10d.o 4_10a1.o 
	$(CC) -g main.o 4_10d.o 4_10a1.o -o 4_10b

4_10c: main.o 4_10c.o 4_10a1.o 4_10d.o
	$(CC) -g main.o 4_10c.o 4_10a1.o -o 4_10c

4_10d.o: 4_10b.c 4_10b0.c
	$(CC) $(CFLAGS) -c -DVER1 4_10b.c
	mv 4_10b.o 4_10d.o

4_10c.o: 4_10b.c 4_10c.c
	$(CC) $(CFLAGS) -c -DVER2 4_10b.c
	mv 4_10b.o 4_10c.o

OUT= 4_10b1.out 4_10b2.out 4_10b3.out 4_10c1.out 4_10c2.out 4_10c3.out
CMP= 4_10b1.cmp 4_10b2.cmp 4_10b3.cmp 4_10c1.cmp 4_10c2.cmp 4_10c3.cmp

4_10b1.out: 4_10b inp1 ; 4_10b < inp1 > 4_10b1.out
4_10b2.out: 4_10b inp2 ; 4_10b < inp2 > 4_10b2.out
4_10b3.out: 4_10b inp3 ; 4_10b < inp3 > 4_10b3.out
4_10c1.out: 4_10c inp1 ; 4_10c < inp1 > 4_10c1.out
4_10c2.out: 4_10c inp2 ; 4_10c < inp2 > 4_10c2.out
4_10c3.out: 4_10c inp3 ; 4_10c < inp3 > 4_10c3.out

test: all $(OUT) $(CMP)
	cmp 4_10b1.out 4_10b1.cmp
	cmp 4_10b2.out 4_10b2.cmp
	cmp 4_10b3.out 4_10b3.cmp
	cmp 4_10c1.out 4_10c1.cmp
	cmp 4_10c2.out 4_10c2.cmp
	cmp 4_10c3.out 4_10c3.cmp
	echo tests done
!EOF!
ls -l makefile
# The following exit is to ensure that extra garbage 
# after the end of the shar file will be ignored.
exit 0