Piotr Krzyżanowski

Poniżej przestawiamy wybrane kody źródłowe programów omawianych w drugiej części podręcznika Obliczenia inżynierskie i naukowe (PWN, 2011). Większość z nich jest dostępna na licencji Creative Commons BY-SA.

Listing 9.1 - harm.c. [Pobierz]

/*
**
** Listing 9.1 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

/* 
Kompilacja: gcc -o harm harm.c
Uruchomienie: ./harm
*/

#include <stdio.h>
#define N 2012

int main(void)
{
	float x;
	unsigned int i;

	x = 0.0;
	for(i = 1; i <= N; i++) 
		x = x + 1.0/i;

	printf("Wartosc sumy x = 1 + 1/2 + .. + 1/%d jest rowna %g\n",  N, x);

	return(0);
}

Listing 9.2 - sinusy.c. [Pobierz]

/*
**
** Listing 9.2 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

/* kompilacja: gcc  -o sinusy sinusy.c -lm */

#include <stdio.h> 
#include <stdlib.h> /* zawiera prototyp funkcji drand48() i srand48() */
#include <math.h>

#define N 15 /* ile liczb wydrukować */

int main(void)
{
int i;
double x, y[N];

srand48(time(NULL)); /* inicjalizacja generatora liczb losowych;
		nie używać, gdy chcemy mieć powtarzalne (sic!) wyniki */

for( i = 0; i < N; i++)
{
	x = 2.0*M_PI*drand48();
	y[i] = sin(x);
	fprintf(stderr, "(%3d) x = %10.5le | sin(x) = %10.5le\n", i, x, y[i]);
}

return(0);
}

Listing 9.3 - zapiszdane.c. [Pobierz]

/*
**
** Listing 9.3 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

/* kompilacja: gcc -o zapiszdane zapiszdane.c -lm */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

const int N = 10; 	/* liczba wierszy */
const int M = 5;	/* liczba kolumn */

#define IJ(i,j,N) (((i)-1) + ((j)-1)*(N))

int main(void)
{
FILE *plik;
double *T;
int i,j;

/* Faza 1: tworzymy tablicę liczb T */

if ((T = (double *) malloc(N*M*sizeof(double))) == NULL)
{
	fprintf(stderr, "Nie mozna przydzielic pamieci\n");
	return(-1);
}

for( j = 1; j <= M; j++)
	for( i = 1; i <= N; i++)
		T[IJ(i,j,N)] = exp(-i/(double)N -j);
		
/* Faza 2: wypisujemy dane z tablicy T do pliku tablica.dat */

if ((plik = fopen("tablica.dat","w")) == NULL)
{
	fprintf(stderr, "Nie mozna otworzyc pliku\n");
	return(-1);
}

fprintf(plik, "%8d %8d\n", N, M);
for( i = 1; i <= N; i++)
{
	for( j = 1; j <= M; j++)
		fprintf(plik, "%16.10E ", T[IJ(i,j,N)]);
	fprintf(plik, "\n");
}

/* Faza 3: sprzątanie */
ferror(plik);
fclose(plik);
free(T);
return(0);
}

Listing 9.4 - wczytajdane.c. [Pobierz]

/*
**
** Listing 9.4 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

/* kompilacja: gcc -o wczytajdane wczytajdane.c */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define IJ(i,j,N) (((i)-1) + ((j)-1)*(N))

int main(int argc, char** argv)
{
FILE *plik;
double *X; /* tu będzie wczytana macierz */
int N,M; /* zmienne na wymiary macierzy */
int i,j;

/* wczytujemy dane z pliku do tablicy - najpierw musimy ustalić jej wymiary */

if ( (plik = fopen(argv[1], "r")) == NULL )
{
	fprintf(stderr, "Nie można otworzyć pliku '%s'\n", argv[1]);
	return(-1);
}

fscanf(plik, "%d %d\n", &N, &M);

/* znając wymiary, alokujemy pamięć na macierz */

if ((X = (double *) malloc(N*M*sizeof(double))) == NULL)
{
	fprintf(stderr, "Nie można przydzielić pamięci\n");
	return(-1);
}

fprintf(stderr, "Wczytuje macierz %d x %d\n", N, M);

for( i = 1; i <= N; i++)
{
	for( j = 1; j <= M; j++)
		fscanf(plik, "%lE", &X[IJ(i,j,N)]);
}

/* ...i od tej pory możemy używać tablicy X, 
  wypełnionej danymi z pliku; ale my już kończymy program... */

fclose(plik);
free(X);
return(0);

}

Listing 9.5 - IJ.h. [Pobierz]

/*
**
** Listing 9.5 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz swobodnie używać tego kodu
**
*/

/* makra dostępu do macierzy rozwiniętej kolumnami */
/* elementy indeksowane od "1" */
#define I(i) ((i)-1) /* wektor */
#define IJ(i,j,n) ((i)-1+((j)-1)*(n)) /* macierz n x m */
#define IJK(i,j,k,n,m) ((i)-1+((j)-1)*(n)+((k)-1)*(n*m)) /* macierz n x m x k */

Listing 9.6 - testblas.c. [Pobierz]

/*
**
** Listing 9.6 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

/* Kompilacja:
	gcc -o testblas testblas.c -lblas -lgfortran -lm 
*/

#include <stdio.h>
#define N 3

extern double dnrm2_(int*,double*,int*); /* prototyp DNRM2 z BLAS */

int main(void)
{
	int n=N, incx=1;
	double x[N]= {3,0,4};
	
	printf("Norma zadanego wektora: %e\n", dnrm2_(&n, x, &incx));
	return(0);
}

Listing 9.7 - mx2.c. [Pobierz]

/*
**
** Listing 9.7 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

#include <stdio.h>
#include <stdlib.h>
#include "mx2.h"

double *mx2alloc(int n, int m)
{
/* przydziela pamięć na macierz n x m liczb typu double */
	return((double *) malloc(n*m*sizeof(double)));
}

void mx2free( double *matrix )
{
/* zwalnia pamięć wskazywaną przez matrix */
	free(matrix);
}

void mx2rand( double *matrix, int n, int m)
{
/* generuje losową macierz n x m o wartościach z przedziału [0,1) */
int i,j;
for( j = 1; j <= m; j++)
	for( i = 1; i <= n; i++)
		matrix[IJ(i,j,n)] = drand48();
}

int mx2write( double *matrix, int n, int m, char *nazwapliku)
{
/* zapisuje (wierszami) macierz do pliku tekstowego;
   w pierwszej linii podajemy wymiary macierzy */
int i,j;
FILE *plik;

if ((plik = fopen(nazwapliku,"w")) == NULL)
{
	fprintf(stderr, "Nie mozna otworzyc pliku: %s\n", nazwapliku);
	return(-1);
}

fprintf(plik, "%8d %8d\n", n, m);
for( i = 1; i <= n; i++)
{
	for( j = 1; j <= m; j++)
		fprintf(plik, "%16.10E ", matrix[IJ(i,j,n)]);
	fprintf(plik, "\n");
}

fclose(plik);
return(0);
}

double mx2frobnorm(double *matrix, int n, int m)
{
/* 
Norma Frobeniusa macierzy dwuwymiarowej: $||A||_F = \sqrt{\sum_{i,j} A_{ij}^2}$ Wystarczy potraktować macierz jako długi wektor (a to już i tak zrobiliśmy), 
a potem obliczyć jego normę euklidesową, korzystając z DNRM2 z biblioteki BLAS 
*/
static int incx=1;
int size;
	size = n*m;
	return(dnrm2_(&size, matrix, &incx)); 
}

Listing 9.8 - mx2.h. [Pobierz]

/*
**
** Listing 9.8 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

/* plik nagłówkowy dla biblioteczki macierzowej: macierze 2-wymiarowe */

#include "IJ.h" /* makro dostępu do elementów macierzy */

/* prototyp funkcji DNRM2 z BLAS */
extern double dnrm2_(int *, double *, int *);

/* prototypy definiowanych funkcji */

double *mx2alloc(int n, int m);
void mx2free( double *matrix );
void mx2rand( double *matrix, int n, int m);
int mx2write( double *matrix, int n, int m, char *nazwapliku);
double mx2frobnorm(double *matrix, int n, int m);

Listing 9.9 - testmx.c. [Pobierz]

/*
**
** Listing 9.9 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

#include <stdio.h>
#include "mx2.h"

int main(void)
{
int N=10, M=5;
double *A;

A = mx2alloc(N,M);
mx2rand(A,N,M);
A[IJ(1,1,N)] = 1.0;
mx2write(A,N,M,"macierzA.dat");
printf("Norma Frobeniusa wygenerowanej macierzy = %e\n", 
	mx2frobnorm(A,N,M));
return(0);
}

Listing 9.10 - Makefile. [Pobierz]

##
## Listing 9.10 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
## (C) Piotr Krzyżanowski, 2011
##
## Możesz swobodnie używać tego kodu
##

# Makefile dla programu testblas.c
CC = gcc
CFLAGS = -O3
CLIBS = -lblas -lgfortran -lm

testblas: testblas.c 
	$(CC) $(CFLAGS) -o $@ $< $(CLIBS)

Listing 9.13 - Makefile-umfpack. [Pobierz]

##
## Listing 9.13 i 10.11 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
## (C) Piotr Krzyżanowski, 2011
##
## Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
## zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
##

# Zapisz ten plik pod nazwą "Makefile"!

#
# "uniwersalny" Makefile przystosowany do kompilacji pliku testumfpacksolve.c
# z listingu 10.10
#

#--------- definicja konkretnego zadania (projektu) ---------
PROJECT = testumfpacksolve

#--------- ogólne definicje ---------
CC = gcc
CFLAGS = -O3 -funroll-loops -static
CLIBS = -lumfpack -lamd -llapack -lblas -lgfortran -lm
CINCDIR = -I/usr/include/suitesparse -I..

#--------- pliki źródłowe i pośrednie -----------	
OBJS = mmread.o mmio.o $(PROJECT).o

#--------- ogólne reguły ---------
%.o: %.c Makefile $(HEADERS)
	$(CC) $(CFLAGS) $(CINCDIR) -c $< -o $@ 

#--------- główne zadanie ---------
$(PROJECT): $(OBJS)
	$(CC) $(CFLAGS) $(CLIBDIR) -o $@ $(OBJS) $(CLIBS) -lrt

#--------- praktyczne zadania ogólne ---------
clean:
	rm $(OBJS) $(PROJECT)

Listing 10.1 - sinusy-gsl.c. [Pobierz]

/*
**
** Listing 10.1 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

/* Kompilacja: 
	gcc -o sinusy-gsl sinusy-gsl.c -lgsl -lgslcblas -lm
*/
#include <stdio.h> 
#include <gsl/gsl_rng.h> /* zawiera prototypy funkcji losowych GSL */
#include <math.h>

#define N 15 /* ile liczb wydrukować */

int main(void)
{
int i;
double x,y;
gsl_rng *rng; /* wskaźnik do generatora, którego będziemy używać */

gsl_rng_env_setup(); 	/* inicjalizacja domyślnego generatora 
			gsl_rng_default i jego ziarna,
			na podstawie zmiennych środowiskowych 
			GSL_RNG_TYPE i GSL_RNG_SEED. 
			Domyślnie GSL_RNG_TYPE=gsl_rng_mt19937
			oraz GSL_RNG_SEED=0 */ 
			
rng = gsl_rng_alloc(gsl_rng_default); /* utworzenie instancji rng 
			generatora liczb losowych domyślnego typu */

gsl_rng_set(rng, time(NULL)); /* zasianie wygenerowanego generatora 
			bieżącą wartością czasu systemowego */

for( i = 0; i < N; i++)
{
	x = 2.0*M_PI*gsl_rng_uniform(rng);
	y = sin(x);
	fprintf(stderr, "(%3d) x = %10.5le sin(x) = %10.5le\n", i, x, y);
}

gsl_rng_free(rng); /* usuwamy niepotrzebny już generator */

return(0);
}

Listing 10.2 - testgslint.c. [Pobierz]

/*
**
** Listing 10.2 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

double F(double X, void * param) /* wrapper dla funkcji sin(x)/x */
{
	return(sin(X)/X);
}


int main(void)
{
	gsl_function f; /* argument z funkcją podcałkową */

	double A,ABSERR,B, EPSABS,EPSREL,RESULT;
	int IER,NEVAL;

	gsl_integration_workspace *workspace;

	int KEY, LIMIT; 

	/* przygotowujemy argument z funkcją podcałkową */
	f.function = &F;
	
	A = 0.0E0; B = 10*M_PI; /* przedział całkowania */
	EPSABS = 0.0E0; EPSREL = 1.0E-3; /* tolerancja błędu */

	/* parametry specyficzne dla QAG */
	KEY = 1; /* tzn. użyj minimalnej liczby punktów kwadratury bazowej */
	LIMIT = 100; /* maksymalny podział przedziału całkowania */
	workspace  = gsl_integration_workspace_alloc(LIMIT);
	
	/* całkujemy: QAG! */

	IER = gsl_integration_qag(&f, A, B, EPSABS, EPSREL, 
					LIMIT, KEY, workspace, &RESULT, &ABSERR);

	if (IER != 0)
		fprintf(stderr,"GSL_QAG: Kłopoty z całkowaniem\n");
	fprintf(stderr,"Całka: %g Est. błąd: %g IER: %d\n", RESULT, ABSERR,  IER);

	gsl_integration_workspace_free(workspace);
	
}

Listing 10.4 - sinusy-acml.c. [Pobierz]

/*
**
** Listing 10.4 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

/* Kompilacja wersji dla x86 (tzn. 32-bitowej): 
	gcc -static -o sinusy-acml sinusy-acml.c \
	 -I $HOME/Software/acml4.3.0/gfortran32/include \
	 -L $HOME/Software/acml4.3.0/gfortran32/lib \
	 -lacml -lgfortran -lm -lrt
Kompilacja wersji dla x86-64 (tzn. 64-bitowej): 
	gcc -static -o sinusy-acml sinusy-acml.c \
	 -D ACML_MV \
	 -I $HOME/Software/acml4.3.0/gfortran64/lib \
	 -L $HOME/Software/acml4.3.0/gfortran64/lib \
	 -lacml_mv -lacml -lgfortran -lm -lrt
*/
#include <stdio.h> 
#include <stdlib.h> 
#include <acml.h> /* zawiera prototypy funkcji losowych ACML */
#ifdef ACML_MV
	#include <acml_mv.h>/* zawiera prototypy funkcji wektorowych ACML */
#endif
#include <math.h>

/* stałe opisujące paramtery generatora Mersenne Twister (MT) w ACML */
#define DRAND_MT_CODE 2 /* kod generatora */
#define DRAND_MT_LSEED 624 /* rozmiar tablicy z ziarnem */
#define DRAND_MT_LSTATE 20 /* rozmiar tablicy stanu */

#define N 15 /* ile liczb wydrukować */

int main(void)
{
int i;
double x[N],y[N];

/* zmienne wymagane do inicjalizacji generatora liczb losowych */
int drand_lseed = DRAND_MT_LSEED; 
int drand_seed[DRAND_MT_LSEED];
int drand_lstate = DRAND_MT_LSTATE;
int drand_state[DRAND_MT_LSTATE];
int drand_info;

/* inicjalizacja ziarna generatora MT */
for(i=0; i < DRAND_MT_LSEED; i++)
	drand_seed[i] = rand();
	
/* inicjalizacja generatora bazowego */
drandinitialize(DRAND_MT_CODE, 1, drand_seed, &drand_lseed, 
		drand_state, &drand_lstate, &drand_info);

/* generujemy od razu N liczb losowych! */
dranduniform(N, 0.0, 2.0*M_PI, drand_state, x, &drand_info); 

#ifdef ACML_MV
	vrda_sin(N, x, y); /* wektorowa instrukcja sin() */
#else
	for( i = 0; i < N; i++)
	{
		y[i] = sin(x[i]);
	}
#endif

for( i = 0; i < N; i++)
	fprintf(stderr, "(%3d) x = %10.5le sin(x) = %10.5le\n", i, x[i], y[i]);

return(0);
}

Listing 10.5 - blaslapack.h. [Pobierz]

/*
**
** Listing 10.5 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

/* przydatne pliki nagłówkowe dla procedur BLAS i LAPACK-a */

#include <stdlib.h>
#include "IJ.h" /* makra dostępu do macierzy i wektorów - zob. listing 9.5 */

/* przydatne stałe */
static int BLASONE = 1, BLASMONE = -1, BLASZERO = 0;

/* BLAS Level 1 */
extern double dnrm2_(int* N, double* X, int* INCX);

/* BLAS Level 2 */
extern void dgemv_( char * TRANS, int * M, int * N, \
	double * ALPHA, double * A,  int * LDA, double * X,  int * INCX, \
	double * BETA, double * Y, int * INCY );

/* BLAS Level 3 */
extern void dgemm_(char * TRANSA,char * TRANSB, int * M, int * N, int * K, \
	double * ALPHA, double * A,  int * LDA,  double * B,  int * LDB, \
	double * BETA, double * C, int * LDC );

/* LAPACK */
extern void dgesv_( int * N, int * NRHS, double * A, int * LDA, \
	int * IPIV, double * B, int * LDB, int * INFO );
extern void dsysv_(char *, int *, int *, double *, int *, \
	int *, double *, int *, double *, int *, int *);
extern void dgels_( char * TRANS,  int * M, int * N, int * NRHS, \
	double * A, int * LDA, double * B, int * LDB, \
	int * WORK, double * LWORK, int * INFO);

Listing 10.6 - testdgemv.c. [Pobierz]

/*
**
** Listing 10.6 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

#include <stdio.h>
#include "blaslapack.h"

int main()
{
	int N = 5, i, j;
	double *A, *x, *y;
	
	/* generujemy przykładową macierz N x N jako odwrotną macierz Laplasjanu */
	A = (double *)malloc(N*N*sizeof(double));
	for (j = 1; j <= N; j++)
	{
		for (i = 1; i <= j; i++)
			A[IJ(i,j,N)] = i;
		for (i = j+1; i <= N; i++)
			A[IJ(i,j,N)] = j;			
	}
	
	/* tworzymy wektory x i y, inicjalizujemy x */
	x = (double *)malloc(N*sizeof(double));
	y = (double *)malloc(N*sizeof(double));
	for (i = 1; i <= N; i++)
		x[I(i)] = (double)i;

	/* obliczamy y = 5*A*x, korzystając z procedury BLAS Level 2: 
	DGEMV dla zadania y = 5*A*x + 0*y */

	{
		char TRANS = 'N'; double ALPHA = 5.0, BETA = 0.0;

		dgemv_(&TRANS, &N, &N, &ALPHA, A,  &N,  x,  &BLASONE, 
			&BETA, y, &BLASONE );
	}
	
	for (i = 1; i <= N; i++)
		printf("%E\n",y[I(i)]);
	
	return(0);
}

Listing 10.10 - testumfpacksolve.c. [Pobierz]

/*
**
** Listing 10.10 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
*/

#include <stdio.h>
#include <malloc.h>
#include "umfpack.h"

#define CHECK(status) 	{if (status != 0) \
	{ fprintf(stderr, "==> An error occured. Aborting program!\n"); \
	 return(-1);}}

/* procedura wczytująca macierz do formatu AIJ indeksowanego od zera */
int mmreadsparse(char *filename, int* N, int* M, int* nz, int** I, int** J,
double** Val );

int main(int argc, char** argv)
{
int M, N, nz; /* wymiary macierzy i liczba niezerowych elementów */  
int *Ap, *Ai; double *Av; /* tu będziemy przechowywać macierz w formacie CSC */
int *I, *J; double *Val; /* tu będziemy przechowywać macierz w formacie AIJ */

double *x, *b; /* wektory prawej strony i rozwiązania */
int i;

/* zmienne używane przez procedury UMFPACK */
double Control[UMFPACK_CONTROL], Info[UMFPACK_INFO];
void *Symbolic, *Numeric;
int status; /* kod zakończenia procedury: różny od zera oznacza kłopoty */

/* ustawiamy parametry sterujące UMFPACKa */
umfpack_di_defaults(Control);    
Control[UMFPACK_PRL] = 2; /* poziom szczegółowości diagnostyki */
umfpack_di_report_control(Control);

if (argc < 2) /* sposób użycia programu */
{
	fprintf(stderr,"Wywołanie: %s PlikMacierzy [PlikPrawejStrony [PlikRozwiązania]]\n", argv[0]);
	return(-1);
}
/* wczytujemy macierz w formacie MatrixMarket do formatu AIJ*/
if (mmreadsparse(argv[1], &N, &M, &nz, &I, &J, &Val ) != 0)
	return(-1);

/* Transformujemy macierz z formatu AIJ do formatu CSC */
Ap = (int *) malloc((N+1) * sizeof(int)); /* N+1, a nie N: to ważne!! */
Ai = (int *) malloc(nz * sizeof(int));
Av = (double *) malloc(nz * sizeof(double));
umfpack_di_triplet_to_col(N,N,nz, I,J,Val, Ap,Ai,Av,(int *)NULL );

umfpack_di_report_matrix(N,N, Ap, Ai, Av, 1, Control);
umfpack_di_report_triplet(N,N,nz, I, J, Val, Control);

/* zwalniamy pamięć zajmowaną przez niepotrzebne zmienne pomocnicze */
free(I); free(J); free(Val);
 
/* faktoryzacja symboliczna, reordering */
status = umfpack_di_symbolic(N, N, Ap, Ai, Av, &Symbolic, Control, Info); 
umfpack_di_report_status(Control,status);
CHECK(status);

/* faktoryzacja numeryczna */
status = umfpack_di_numeric(Ap, Ai, Av, Symbolic, &Numeric, Control, Info); 
umfpack_di_report_status(Control,status);
CHECK(status);

/* szykujemy wektory rozwiązania i prawej strony */
x = (double *)calloc(N,sizeof(double));
b = (double *)calloc(N,sizeof(double));  /* calloc zeruje alokowaną pamięć */

/* wczytujemy wektor prawej strony (o ile jest podany) w formacie MatrixMarket*/
if (argc > 2) 
{
	int Nrhs, Mrhs;
	
	if (mmreadsparse(argv[2], &Nrhs, &Mrhs, &nz, &I, &J, &Val) != 0)
		return(-1);
	else
		fprintf(stderr,"Wczytano prawą stronę : %dx%d\n", Nrhs, Mrhs);
	if (Mrhs > 1)
	{
		fprintf(stderr,"Oczekiwano jednego wektora prawej strony!\n");
		return(-1);
	}
	
	/* zapisujemy go do wektora b */
	for (i = 1; i <= nz; i++)
		b[I[i]] += Val[i];

	/* zwalniamy pamięć zajmowaną przez niepotrzebne zmienne pomocnicze */
	free(I); free(J); free(Val);
}

/* rozwiązanie układu Ax=b */
status = umfpack_di_solve(UMFPACK_A, Ap, Ai, Av, x, b, Numeric, Control, Info); 
umfpack_di_report_status(Control,status);
CHECK(status);

/* statystyki */
umfpack_di_report_info(Control,Info);

/* sprzątanie */
umfpack_di_free_symbolic(&Symbolic);
umfpack_di_free_numeric(&Numeric);

/* ... korzystamy z wyznaczonego rozwiązania x .... */
/* zapisujemy je do pliku (tekstowego), o ile było takie życzenie */
if (argc > 3) 
{
	FILE *plik;
	if ((plik = fopen( argv[3], "w")) == NULL)
	{
		fprintf(stderr, "Nie mozna zapisac pliku %s\n", argv[3]);
		return(-1);
	}
	for (i = 0; i < N; i++)
		fprintf(plik, "%26.10e \n", x[i]);
	fclose(plik);	
}

/* zwalniamy pamięć zajmowaną przez niepotrzebne zmienne */
free(Ap); free(Ai); free(Av); free(x); free(b);

return(0);
}

Listing 11.1 - timer.c. [Pobierz]

/*
**
** Listing 11.1 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
** (C) Piotr Krzyżanowski, 2011
**
** Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
** zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
**
*/

/* Stoper w C, mierzący czas na trzy różne sposoby, 
za pomocą funkcji clock(), times(), clock_gettime(). 
Sposób użycia:
	tic();
	... instrukcje ...
	czas = toc();
*/

#include <stdio.h>
#include <unistd.h>	/* zawiera pewne stałe systemowe */
#include <sys/times.h> 	/* zawiera funkcję times() */
#include <time.h> 	/* zawiera funkcję clock_gettime() i clock() */

static clock_t start_clk, stop_clk; /* dla clock() */
static struct tms start_tms, stop_tms; /* dla times() */
static struct timespec start_tmsp, stop_tmsp; /* dla clock_gettime() */

double tic(void)
{
	fprintf(stderr, ">>>> Stoper START\n");
	
	start_clk = clock();
	times(&start_tms);
	clock_gettime(CLOCK_REALTIME, &start_tmsp);
	
	return(0.0);
}

double toc(void)
{
double cpu_timer_clk;
double cpu_timer_total;
double cpu_timer_real;

	stop_clk = clock();
	times(&stop_tms);
	clock_gettime(CLOCK_REALTIME, &stop_tmsp);
	
	cpu_timer_total = (stop_tms.tms_utime-start_tms.tms_utime) *
			(1.0/sysconf(_SC_CLK_TCK));
	cpu_timer_clk = (stop_clk-start_clk) *
			(1.0/CLOCKS_PER_SEC);
	cpu_timer_real = stop_tmsp.tv_sec-start_tmsp.tv_sec +
			1e-9*(stop_tmsp.tv_nsec-start_tmsp.tv_nsec);
	
	fprintf(stderr, "<<<< Stoper STOP: %g (TOTAL CPU) %g (CLOCK) %g (REAL) sekund\n", 
			cpu_timer_total, cpu_timer_clk, cpu_timer_real);
	
	return(cpu_timer_real);
}