PDF de programación - Implementación en fortran 77 del método de ortonormalización de Gram–Schmidt en R

Implementación en fortran 77 del método de

ortonormalización de Gram–Schmidt en Rn

Pablo Santamaría

v0.2 (Marzo 2008)

Sea S el subespacio de Rn generado por el conjunto de m vectores linealmente independientes
{v1, v2, . . . , vm} (m ≤ n). El método de ortonormalización de Gram-Schmidt construye una base
ortornomal {w1, w2, . . . , wm} para S a partir de la base anterior según el procedimiento:

w1 =

wj =

,

v1
kv1k
vj − Pj−1

vj − Pj−1


l=1 hwl|vj iwl
l=1 hwl|vj iwl


j = 2, . . . , m,

siendo hv|ui = vtu el producto interno canónico de Rn y kvk = phv|vi la norma asociada. Una
implementación algorítmica de este procedimiento puede plantearse como sigue.

Dados v1, v2, . . . , vm
Para j = 1, 2, . . . , m

Tomar wj = vj
Para l = 1, 2 . . . , j − 1

Calcular rlj = hwl|vji
Tomar wj = wj − rlj wl

Calcular rjj = kwjk
Tomar wj =

wj
rjj

Sin embargo, esta implementación directa del método no constituye una buena implementación
numérica, por cuanto al utilizar aritmética de punto flotante los vectores sucesivos calculados se
apartan de la ortogonalidad debido a los errores de redondeo involucrados. Una implementación
matemáticamente equivalente, pero numéricamente más estable procede como sigue

Dados v1, v2, . . . , vm
Para j = 1, 2, . . . , m

Tomar wj = vj
Para l = 1, 2 . . . , j − 1

Calcular rlj = hwl|wji
Tomar wj = wj − rlj wl

Calcular rjj = kwjk
Tomar wj =

wj
rjj

Para la implementación computacional de este procedimiento disponemos los m vectores vj,
cada uno de n componentes, como las columnas de una matriz A cuyas dimensiones lógicas son,
pues, n×m. Como es usual, la dimensión física de la matriz será nmax×nmax. La siguiente subrutina
implementa el método. Obsérvese que los vectores ortonormalizados vuelven en la misma matriz
A, cuyos valores originales, son pues, destruidos.

1

SUBROUTINE GRAMSCHMIDT(A,N,M,NMAX)

IMPLICIT NONE
INTEGER NMAX,N,M
DOUBLE PRECISION A(NMAX,*)
INTEGER I,J,L
DOUBLE PRECISION PI

DO J=1,M

DO L=1,J-1

----------------------------
CALCULAR PI = <W_L|W_J>
----------------------------
PI = 0.0D0
DO I=1,N

PI = PI + A(I,L)*A(I,J)

END DO
----------------------------
TOMAR W_J = W_J - PI*W_J
----------------------------
DO I=1,N

A(I,J) = A(I,J)-PI*A(I,L)

END DO

END DO
----------------------------
CALCULAR NORMA DE W_J
----------------------------
PI = 0.0D0
DO I=1,N

PI = PI + A(I,J)**2

END DO
PI = SQRT(PI)
----------------------------
NORMALIZAR W_J
----------------------------
DO I=1,N

A(I,J) = A(I,J)/PI

END DO

END DO
RETURN
END

*
*
*

*
*
*

*
*
*

*
*
*

Espero que este apunte les haya resultado útil, y como siempre, cualquier duda o sugerencia

envíenla por e-mail a [email protected].

2