Example 1. Gauss elimination algorithm
PROGRAM GAUSS C Solving linear equation system Ax = b PARAMETER ( N = 100 ) REAL A( N, N+1 ), X( N ) C A : Coefficient matrix with dimension (N,N+1) C Right hand side vector of linear equations is stored C in last (N+1)-th column of matrix A C X : Vector of unknowns C N : Number of linear equations *HPF$ DISTRIBUTE A (BLOCK,*) *HPF$ ALIGN X(I) WITH A(I,N+1) C Initialization *HPF$ INDEPENDENT DO 100 I = 1, N DO 100 J = 1, N+1 IF (( I .EQ. J ) THEN A(I,J) = 2.0 ELSE IF ( J .EQ. N+1) THEN A(I,J) = 0.0 ENDIF ENDIF 100 CONTINUE C C Elimination C DO 1 I = 1, N *HPF$ INDEPENDENT DO 5 J = I+1, N DO 5 K = I+1, N+1 A(J,K) = A(J,K) - A(J,I) * A(I,K) / A(I,I) 5 CONTINUE 1 CONTINUE C First calculate X(N) X(N) = A(N,N+1) / A(N,N) C C Solve X(N-1), X(N-2), ...,X(1) by backward substitution C DO 6 J = N-1, 1, -1 *HPF$ INDEPENDENT DO 7 I = 1, J A(I,N+1) = A(I,N+1) - A(I,J+1) * X(J+1) 7 CONTINUE X(J) = A(J,N+1) / A(J,J) 6 CONTINUE PRINT *, X END
PROGRAM JACOBI PARAMETER (K=8, ITMAX=20) REAL A(K,K), B(K,K) *HPF$ DISTRIBUTE A (BLOCK, BLOCK) *HPF$ ALIGN B(I,J) WITH A(I,J) C arrays A and B with block distribution PRINT *, '******** TEST_JACOBI_HPF ********' C nest of two independent loops, iteration (i,j) will be executed C on processor, which is owner of element A(i,j) *HPF$ INDEPENDENT DO 1 J = 1, K *HPF$ INDEPENDENT DO 1 I = 1, K A(I,J) = 0. IF(I.EQ.1 .OR. J.EQ.1 .OR. I.EQ.K .OR. J.EQ.K) THEN B(I,J) = 0. ELSE B(I,J) = 1. + I + J ENDIF 1 CONTINUE DO 2 IT = 1, ITMAX *HPF$ INDEPENDENT DO 21 J = 2, K-1 *HPF$ INDEPENDENT DO 21 I = 2, K-1 A(I,J) = B(I,J) 21 CONTINUE *HPF$ INDEPENDENT DO 22 J = 2, K-1 *HPF$ INDEPENDENT DO 22 I = 2, K-1 B(I,J) = (A(I-1,J) + A(I,J-1) + A(I+1,J) + A(I,J+1)) / 4 22 CONTINUE 2 CONTINUE 3 OPEN (3, FILE='JACH.DAT', FORM='FORMATTED') WRITE (3,*) B CLOSE (3) END
Example 3. Red-black successive over-relaxation
PROGRAM REDBLACK
PARAMETER (N1 = 20,N2 = 10)
REAL A(N1,N2),W
INTEGER ITMAX
*HPF$ DISTRIBUTE (BLOCK,BLOCK) :: A
ITMAX = 20
W = 0.5
*HPF$ INDEPENDENT
DO 1 J = 1,N2
*HPF$ INDEPENDENT
DO 1 I = 1,N1
IF (I.EQ.J) THEN
A(I,J) = N1+2
ELSE
A(I,J) = (-(1.))
ENDIF
1 CONTINUE
DO 2 IT = 1,ITMAX
*HPF$ INDEPENDENT
DO 21 J = 1,N2/2-1
*HPF$ INDEPENDENT
DO 21 I = 1,N1/2-1
A(2*I+1,2*J+1) = W/4*(A(2*I,2*J+1)+A(2*I+2,2*J+1)+
+ A(2*I+1,2*J)+A(2*I+1,2*J+2))+(1-W)*A(2*I+1,2*J+1)
21 CONTINUE
*HPF$ INDEPENDENT DO 22 J = 1, N2/2-1
*HPF$ INDEPENDENT
DO 22 I = 1,N1/2-1
A(2*I,2*J) = W/4*(A(2*I-1,2*J)+A(2*I+1,2*J)+A(2*I,2*J-1)+
+ A(2*I,2*J+1))+(1-W)*A(2*I,2*J)
22 CONTINUE
*HPF$ INDEPENDENT
DO 23 J = 1,N2/2-1
*HPF$ INDEPENDENT
DO 23 I = 1,N1/2-1
A(2*I,2*J+1) = W/4*(A(2*I-1,2*J+1)+A(2*I+1,2*J+1)+A(2*I,2*J)
+ +A(2*I,2*J+2))+(1-W)*A(2*I,2*J+1)
23 CONTINUE
*HPF$ INDEPENDENT
DO 24 J = 1,N2/2-1
*HPF$ INDEPENDENT
DO 24 I = 1,N1/2-1
A(2*I+1,2*J) = W/4*(A(2*I,2*J)+A(2*I+2,2*J)+A(2*I+1,2*J-1)+
+ A(2*I+1,2*J+1))+(1-W)*A(2*I+1,2*J)
24 CONTINUE
PRINT *,'IT= ',IT
2 CONTINUE
END
Example 4. Gauss Elimination with Pivoting and Error Correction
***GAUSS ELIMINATION METHOD************************************************
* THE PROGRAM FOR SYSTEM OF LINEAR EQUATION BY GAUSS ELIMINATION METHOD *
* WITH COMPLETE PIVOTING AND ERROR CORRECTION *
***************************************************************************
*
program GAUSS
integer nd,i,j,n,m,nn,k,kk,p,r,index,determ,chec,z
parameter (nd = 40)
doubleprecision AM(nd,nd),XM(nd),CM(nd),EM(nd),PE(nd),DX(nd),
$XX(nd),DM(nd),ORDC(nd),YM(nd),qt,d,t
parameter (epsil = 1.0e-9)
**
* SEGMENT FOR READING THE EQUATION
**
write(6,*) ' Enter number of equations, n:'
read(*,*) n
m = n + 1
write(6,*) ' Enter augmented coefficient matrix row-wise:'
read(*,*) ((AM(i,j),j = 1,m),i = 1,n)
**
z = 0
nn = n - 1
chec = 1
**
* ESTABLISH INITIAL ORDER IN THE COLUMN ORDER VECTOR
**
do 270 i = 1,n
ORDC(i) = i
270 continue
**
* SEGMENT FOR PARTIAL PIVOTING
**
do 300 p = 1,nn
CALL PIVOT(AM,n,ORDC,p,ORD)
**
* TRIANGULARIZATION BY ELIMINATING THE VARIABLES
**
kk = p + 1
do 290 i = kk,n
if(abs(AM(p,p)).lt.epsil) then
chec = 0
else
qt = AM(i,p)/AM(p,p)
endif
do 280 j = p,m
AM(i,j) = AM(i,j) - qt * AM(p,j)
280 continue
290 continue
300 continue
**
* CHECKING FOR THE SINGULARITY OF COEFFICIENT MATRIX
**
determ = 1
do 310 i = 1,n
if ((abs(AM(i,i)).lt.epsil).or.(chec.eq.0)) then
determ = 0
write(6,*)'--------------------------------------------------'
write(6,*)'Coefficient matrix is singular. No Solution exist.'
write(6,*)'--------------------------------------------------'
endif
310 continue
if (determ.eq.1.and.chec.eq.1) then
**
* BACK SUBSTITUTION
**
XM(n) = AM(n,m)/AM(n,n)
do 330 i = nn,1,-1
sum = 0.0
k = i + 1
do 320 j = k,n
sum = sum + AM(i,j) * XM(j)
320 continue
XM(i) = (AM(i,m) - sum)/AM(i,i)
330 continue
**
* REARRANGE THE SOLUTION VECTOR
**
do 340 i = 1,n
j = ORDC(i)
YM(j) = XM(i)
340 continue
**
* OUTPUTTING THE SOLUTION VECTOR
**
write(6,*)'--------------------------------------------------'
write(6,*) ' The solution Vector is:'
write(6,*)
write(6,341) (YM(i),i = 1,n)
341 format('',30x,f12.4)
write(6,*)'--------------------------------------------------'
**
* CHECK FOR THE VALUE OF THE DETERMINANT
**
d = 1.0
do 350 i = 1,n
d = d * AM(i,i)
350 continue
d = d * (-1) **z
write(6,351) 'The determinant of coefficient Matrix is:',d
351 format('0',a41,2x,f10.3)
write(6,*)'--------------------------------------------------'
write(6,*) ' ENTER [RETURN] TO CONTINUE:'
read(*,*)
**
* SUBSTITUTION OF SOLUTION IN THE EQUATION TO GET THE RIGHT HAND SIDE
**
do 370 i = 1,n
CM(i) = 0.0
do 360 index = 1,n
CM(i) = CM(i) + AM(i,index) * XM(index)
360 continue
370 continue
**
* CALCULATION OF RIGHT HAND SIDE ERROR VECTOR
**
do 380 i = 1,n
EM(i) = AM(i,m) - CM(i)
380 continue
write(6,*)'---------------------------------------------------'
write(6,*) ' The right-hand side error vector is:'
write(6,381) (EM(i),i = 1,n)
381 format('',28x,f12.4)
write(6,*)'---------------------------------------------------'
write(6,*) 'The relative percentage errors of the right-hand sid
$e contants are:'
**
* CALCULATION OF PERCENTAGE ERROR OF RIGHT HAND SIDE CONSTANT
**
do 390 i = 1,n
PE(i) = EM(i)/CM(i)
390 continue
CALL pr(PE,n,1)
write(6,*)
**
**
do 400 i = 1,n
AM(i,m) = EM(i)
400 continue
**
* DETERMINATION OF CORRECTION FACTOR VECTOR
**
DX(n) = AM(n,m)/AM(n,n)
do 420 i = nn,1,-1
sum = 0.0
k = i + 1
do 410 j = k,n
sum = sum + AM(i,j) * DX(j)
410 continue
DX(i) = (AM(i,m) - sum)/AM(i,i)
420 continue
write(6,*)'---------------------------------------------------'
write(6,*) ' The correction factor vector for variable is:'
CALL pr(DX,n,1)
write(6,*)
**
write(6,*)'---------------------------------------------------'
write(6,*) ' The corrected solution vector is:'
write(6,*)
do 430 i = 1,n
XX(i) = XM(i) - DX(i)
430 continue
**
* REARRANGE THE SOLUTION VECTOR
**
do 440 i = 1,n
j = ORDC(i)
DM(j) = XX(i)
440 continue
write(6,441) (DM(i),i = 1,n)
441 format('',28x,f12.4)
write(6,*)'---------------------------------------------------'
endif
stop
end
**********************************************************************
*
***COMPLETE PIVOTING (Loop in the calling program)********************
* *
subroutine PIVOT(A,n,ORDC,i,ORD)
integer nd,n,m,i,j,ii,jj,col,row,p,r
parameter (nd = 40)
doubleprecision A(nd,nd),ORD(nd),ORDC(nd),t,tem
**
* COMPLETE PIVOTING
**
row = i
col = i
do 30 p = i,n
do 20 r = i,n
if (abs(A(row,col)).lt.abs(A(p,r))) then
row = p
col = r
endif
20 continue
30 continue
**
if (col.ne.i) then
do 40 ii = 1,n
tem = A(ii,i)
A(ii,i) = A(ii,col)
A(ii,col) = tem
40 continue
t = ORDC(i)
ORDC(i) = ORDC(col)
ORDC(col) = t
endif
*
m = n + 1
if (row.ne.i) then
do 50 jj = 1,m
tem = A(i,jj)
A(i,jj) = A(row,jj)
A(row,jj) = tem
50 continue
t = ORD(i)
ORD(i) = ORD(row)
ORD(row) = t
endif
*
return
end
***SUBROUTINE FOR PRINTING THE MATRIX*********************************
* *
subroutine pr(C,m,n)
integer nd,m,n,i,j
parameter (nd = 40)
doubleprecision C(nd,nd)
if(n.eq.1) then
do 20 i = 1,m
write(6,10)C(i,1)
10 format('',f10.4)
20 continue
elseif(n.eq.2) then
do 40 i = 1,m
write(6,30)C(i,1),C(i,2)
30 format('',2f10.4)
40 continue
elseif(n.eq.3) then
do 60 i = 1,m
write(6,50)C(i,1),C(i,2),C(i,3)
50 format('',3f10.4)
60 continue
elseif(n.eq.4) then
do 80 i = 1,m
write(6,70)C(i,1),C(i,2),C(i,3),C(i,4)
70 format('',4f10.4)
80 continue
elseif(n.eq.5) then
do 100 i = 1,m
write(6,90)C(i,1),C(i,2),C(i,3),C(i,4),C(i,5)
90 format('',5f10.4)
100 continue
elseif(n.eq.6) then
do 120 i = 1,m
write(6,110)C(i,1),C(i,2),C(i,3),C(i,4),C(i,5),C(i,6)
110 format('',6f10.4)
120 continue
elseif(n.eq.7) then
do 140 i = 1,m
write(6,130)C(i,1),C(i,2),C(i,3),C(i,4),C(i,5),C(i,6),C(i,7)
130 format('',7f10.4)
140 continue
elseif(n.eq.8) then
do 160 i = 1,m
write(6,150)C(i,1),C(i,2),C(i,3),C(i,4),C(i,5),C(i,6),C(i,7),
$C(i,8)
150 format('',8f10.4)
160 continue
elseif(n.eq.9) then
do 180 i = 1,m
write(6,170)C(i,1),C(i,2),C(i,3),C(i,4),C(i,5),C(i,6),C(i,7),
$C(i,8),C(i,9)
170 format('',9f10.4)
180 continue
elseif(n.eq.10) then
do 200 i = 1,m
write(6,190)C(i,1),C(i,2),C(i,3),C(i,4),C(i,5),C(i,6),C(i,7),
$C(i,8),C(i,9),C(i,10)
190 format('',10f10.4)
200 continue
endif
return
end
Example 5. Runge Kutta 4 Method
PROGRAM D15R2
C Driver for routine RKDUMB
PARAMETER(NVAR=2)
DIMENSION VSTART(NVAR)
COMMON /PATH/X(200),Y(10,200)
OPEN(6,FILE='RK.OUT',FORM='FORMATTED')
WRITE(*,*) ' X1=? X2=? NSTEP=?'
READ(*,*) X1,X2,NSTEP
VSTART(1)=1.
VSTART(2)=0.
CALL RKDUMB(VSTART,X1,X2,NSTEP)
WRITE(*,'(/1X,T9,A,T17,A,T31,A/)') 'X','X1 VALUE','X2 VALUE'
WRITE(6,'(/1X,T9,A,T17,A,T31,A/)') 'X','X1 VALUE','X2 VALUE'
DO 11 I=1,NSTEP+1
WRITE(*,'(1X,F10.4,2X,2F12.6)') X(i),Y(1,i),Y(2,i)
WRITE(6,'(1X,F10.4,2X,2F12.6)') X(i),Y(1,i),Y(2,i)
11 CONTINUE
END
C
SUBROUTINE DERIVS(X,Y,DYDX)
PARAMETER(NVAR=2)
DIMENSION Y(NVAR),DYDX(NVAR)
DYDX(1)=Y(2)
DYDX(2)=-Y(2)*x-y(1)
RETURN
END
C
SUBROUTINE RKDUMB(VSTART,X1,X2,NSTEP)
PARAMETER (NMAX=10,NVAR=2)
COMMON /PATH/ XX(200),Y(10,200)
DIMENSION VSTART(NVAR),V(NMAX),DV(NMAX)
DO 11 I=1,NVAR
V(I)=VSTART(I)
Y(I,1)=V(I)
11 CONTINUE
XX(1)=X1
X=X1
H=(X2-X1)/NSTEP
DO 13 K=1,NSTEP
CALL DERIVS(X,V,DV)
CALL RK4(V,DV,X,H,V)
IF(X+H.EQ.X)PAUSE 'Stepsize not significant in RKDUMB.'
X=X+H
XX(K+1)=X
DO 12 I=1,NVAR
Y(I,K+1)=V(I)
12 CONTINUE
13 CONTINUE
RETURN
END
C
SUBROUTINE RK4(Y,DYDX,X,H,YOUT)
PARAMETER (NMAX=10,NVAR=2)
DIMENSION Y(NVAR),DYDX(NVAR),YOUT(NVAR),YT(NMAX),DYT(NMAX),
@DYM(NMAX)
HH=H*0.5
H6=H/6.
XH=X+HH
DO 11 I=1,NVAR
YT(I)=Y(I)+HH*DYDX(I)
11 CONTINUE
CALL DERIVS(XH,YT,DYT)
DO 12 I=1,NVAR
YT(I)=Y(I)+HH*DYT(I)
12 CONTINUE
CALL DERIVS(XH,YT,DYM)
DO 13 I=1,NVAR
YT(I)=Y(I)+H*DYM(I)
DYM(I)=DYT(I)+DYM(I)
13 CONTINUE
CALL DERIVS(X+H,YT,DYT)
DO 14 I=1,NVAR
YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.*DYM(I))
14 CONTINUE
RETURN
END
Example 6. Eigenvalues and Eigenvectors
c Main program for finding the eigenvalues and eigenvectors
c of a symmetric matrix AO (NxN) with Jacobi's method. ITMAX
c and EPS are the stopping criteria, namely the maximum number
c of iterations and error tolerance,respectively.
DIMENSION AO(3,3), A(3,3)
DATA N,ITMAX, EPS /3,100,1.0E-06/
AO(1,1)=2.0
AO(1,2)=8.0
AO(1,3)=10.
AO(2,1)=8.0
AO(2,2)=3.0
AO(2,3)=4.0
AO(3,1)=10.
AO(3,2)=4.
AO(3,3)=7.0
CALL JACOBI(AO,N,A,EPS,ITMAX)
PRINT 10, (AO(I,I),I=1,N),((A(I,J), J=1,N), I=1,N)
10 FORMAT (1X,'EIGENVALUES ARE',/,3X,3E18.8,//3X,
# 'EIGENVECTORS',/,8X,' FIRST',15X,' SECOND',13X,'THIRD',
# /,(3X,3E18.8))
STOP
END
C **********************************
SUBROUTINE JACOBI (AO,N,A,EPS,ITMAX)
DIMENSION AO(N,N), A(N,N)
ITER=0
DO 10 I=1,N
DO 10 J=1,N
A(I,J)=0.0
10 A(I,I)=1.0
20 Z=0.
NM1=N-1
DO 30 I=1,NM1
IP1=I+1
DO 30 J=IP1,N
IF (ABS(AO(I,J)) .LE. Z) GO TO 30
Z=ABS(AO(I,J))
IROW=I
ICOL=J
30 CONTINUE
IF (ITER .EQ. 0) Y=Z*EPS
IF (Z .LT. Y) GO TO 200
DIF=AO(IROW,IROW)-AO(ICOL,ICOL)
TANG= (-DIF+SQRT(DIF**2+4.0*Z**2))/(2.0*AO(IROW,ICOL))
COSE=1.0/SQRT(1.0+TANG**2)
SINE=COSE*TANG
DO 40 I=1,N
ZZ=A(I,IROW)
A(I,IROW)=COSE*ZZ+SINE*A(I,ICOL)
40 A(I,ICOL)=COSE*A(I,ICOL)-SINE*ZZ
I=1
50 IF (I .EQ. IROW) GO TO 60
YY=AO(I,IROW)
AO(I,IROW)=COSE*YY+SINE*AO(I,ICOL)
AO(I,ICOL)=COSE*AO(I,ICOL)-SINE*YY
I=I+1
GO TO 50
60 I=IROW+1
70 IF (I .EQ. ICOL) GO TO 80
YY=AO(IROW,I)
AO(IROW,I)=COSE*YY+SINE*AO(I,ICOL)
AO(I,ICOL)=COSE*AO(I,ICOL)-SINE*YY
I=I+1
GO TO 70
80 I=ICOL+1
90 IF (I .GT. N) GO TO 100
ZZ=AO(IROW,I)
AO(IROW,I)=COSE*ZZ+SINE*AO(ICOL,I)
AO(ICOL,I)=COSE*AO(ICOL,I)-SINE*ZZ
I=I+1
GO TO 90
100 YY=AO(IROW,IROW)
AO(IROW,IROW)=YY*COSE**2+AO(IROW,ICOL)*2.0*COSE*SINE+
# AO(ICOL,ICOL)*SINE**2
AO(ICOL,ICOL)=AO(ICOL,ICOL)*COSE**2+YY*SINE**2-AO(IROW,ICOL)
# *2.0*COSE*SINE
AO(IROW,ICOL)=0.0
ITER=ITER+1
IF (ITER .LT. ITMAX) GO TO 20
200 RETURN
END
Example 7 Matrix Multiplication :
#include<stdio.h>
#include<conio.h>
void main()
{
int mat1[2][2],mat2[2][2],mat3[2][2],i,j,k;
printf("************A program to find the multiplication of the given matrices***********\n\n");
printf("************Enter the value for the 1st matrix************\n");
for(i=0;i<2;i++)
{
for(j=0;j<2;j++)
{
printf("\nEnter the value for %d row,%d column:",i+1,j+1);
scanf("%d",&mat1[i][j]);
}
}
printf("\n************Enter the value for the 2nd matrix************\n");
for(i=0;i<2;i++)
{
for(j=0;j<2;j++)
{
printf("\nEnter the value for %d row,%d column:",i+1,j+1);
scanf("%d",&mat2[i][j]);
}
}
printf("\n****************************The given matrices are*****************************\n\n");
for(i=0;i<2;i++)
{
for(j=0;j<2;j++)
{
printf("%d\t",mat1[i][j]);
}
printf("\n");
}
printf("\n\n");
for(i=0;i<2;i++)
{
for(j=0;j<2;j++)
{
printf("%d\t",mat2[i][j]);
}
printf("\n");
}
printf("\n****************************The multiplication of the matrices is*************\n\n");
for(i=0;i<2;i++)
{
for(j=0;j<2;j++)
{
mat3[i][j]=0;
for(k=0;k<2;k++)
{
mat3[i][j]=mat3[i][j]+(mat1[i][k]*mat2[k][j]);
}
}
}
for(i=0;i<2;i++)
{
for(j=0;j<2;j++)
{
printf("%d\t",mat3[i][j]);
}
printf("\n");
}
getch();
getch();
}
Example 8 Finding Roots(Real and complex) of a Third Order Equation :
#include<conio.h>
#include<stdio.h>
#include<math.h>
void main(void)
{
clrscr();
float a,b,c,d,i,j,k,x1,x2,l,r,z,fxl,u,fxr;
float xl = -90.0;
float xu = 100.0;
float xr;
printf("Enter a, b, c & d:\n");
scanf("%f %f %f %f", &a, &b, &c,&d);
for(i=0;i<1000;i++)
{
fxl = (a*pow(xl,3)) + (b*pow(xl,2) ) + (c*xl) + d;
xr = (xl + xu)/2;
fxr = (a*pow(xr,3)) + (b*pow(xr,2) ) + (c*xr) + d;
if((fxl*fxr)<0)
xu = xr;
else
xl = xr;
if(fxr<0.000001 && fxr>-0.000001)
{
printf("The three roots are as follows: -\n");
printf("x1 = %.2f\n", xr);
break;
}
}
i = b + (xr*a);
j = c + (xr*i);
k = i*i - (4*a*j);
if( k >= 0)
{
x1 = (-i + sqrt(k))/(2*a);
x2 = (-i - sqrt(k))/(2*a);
printf("x2 = %.2f\n",x1);
printf("x3 = %.2f",x2);
}
else
{
l = -k;
r = -i/(2*a);
z = (sqrt(l))/(2*a);
printf("x2 = %.2f + j%.2f\n", r,z);
printf("x3 = %.2f - j%.2f\n", r,z);
}
getch();
}