1! Gilles GRASSEAU (CNRS/IDRIS - France) <Gilles.Grasseau@idris.fr> - Nov. 1996 2 3program test_psgetrs 4!------------------------------------------------------------------- 5! Description : 6! Exemple de factorisation LU puis de resolution du systeme : 7! A * x = b , ou A est une matrice (5,5). 8! 9! Principe : 10! Les matrice/vecteurs A, x et b sont globaux et les matrices/ 11! vecteurs dist_a, dist_x et dist_b sont respectivement les 12! matrices/vecteurs distribuees sur la grille de proc. 2x2. 13! La repartition de la matrice se fait par blocs cycliques 2x2 14! (2 elements par ligne et 2 elements par colonne - 15! voir Exemple de repartition d'une matrice). 16! La repartition des vecteurs se fait par blocs cycliques 2x1 17! (le processeur logique (0,0) possede les blocs {x(1:2) ,x(5)} 18! et le processeur logique (1,0) possede le bloc x(3:4) ). 19! 20! Algorithme : 21! 1 - Initialisation du BLACS et autres. 22! 2 - Distribution de la matrice A et du vecteur x vers les blocs 23! locaux dist_a et dist_x de chaque processeur logique. 24! (voir Distribution de matrice) 25! On a simule la situation ou seul le processeur (0,0) possede 26! la matrice et les vecteurs globaux A, x et b. 27! 3 - Factorisation LU et resolution. 28! 4 - Recuperation et impression des resultats par le processeur 29! logique (0,0). 30! 31! Remarque : 32! 1 - Le tableau de travail 'work' a ete dimensionne a la valeur 33! donnee par INITBUFF pour psgetrs ( > au buffer de psgetrf). 34! 2 - La liste des pivots locaux 'dist_piv' est dimensionnee a 35! 'n_max' + 'rb_size'. 36! 37!------------------------------------------------------------------- 38 implicit none 39 40 integer :: pe, npes ! Identificateur du processeur et 41 ! nombre de processeurs physiques. 42 43 integer, parameter:: nprow=2, npcol=2 ! Grille de proc. logiques. 44 integer :: prow, pcol ! Coord. de mon proc. logique 45 46 integer :: icntxt ! Contexte (grille de procs) 47 integer, dimension(8) :: desc_a ! Descripteur de la matrice. 48 integer, dimension(8) :: desc_x ! Descripteur des vecteurs. 49 50 integer, parameter :: n=5 ! Ordre matrices/vecteurs. 51 real, dimension(n,n) :: a ! Matrice globale a. 52 real, dimension(n) :: x ! Vecteur resultat x(global) 53 real, dimension(n) :: b ! Vecteur b (global) 54 55 integer, parameter :: n_max=3 ! Dim. des matrices/vecteurs 56 ! locaux. 57 real, dimension(n_max,n_max):: dist_a ! Matrice locale. 58 real, dimension(n_max) :: dist_x ! Vecteur local. 59 real, dimension(n_max) :: dist_b ! Vecteur local. 60 61 integer, parameter :: rb_size=2, & ! Taille des blocs lignes 62 cb_size=2 ! et colonnes. 63 64 integer, dimension(n_max+rb_size):: dist_piv ! Vecteur local des 65 ! pivots. 66 67 integer :: info, i, j 68 69 data a/ 0.5, 0.0, 2.3, 0.0,-2.6, & ! Coefficients de la matrice 70 0.0, 0.5,-1.4, 0.0,-0.7, & ! globale A. 71 2.3,-1.4, 0.5, 0.0, 0.0, & 72 0.0, 0.0, 0.0, 1.0, 0.0, & 73 -2.6,-0.7, 0.0, 0.0, 0.5/ 74 75 data b/-2.1, 0.3, 1.5, 0.0,-1.5/ ! Coefficients du vecteur global 76 77 78 integer,parameter:: minb=rb_size, maxb=minb ! Calcul de l'espace 79 integer,parameter:: minp=nprow, maxd=n ! de travail 80 integer,parameter:: wsize=2*maxb*(((maxd/minb)/minp)*maxb+maxb)*8 81 real, dimension(wsize/8) :: work 82 83! Initialisation BLACS et autres. 84!-------------------------------------- 85 call initbuff(work,wsize) ! Init. des buffers internes au PBLAS 86 call blacs_pinfo(pe, npes) ! Init. mon PE, nbre de procs physiques 87 call blacs_gridinit(icntxt, & ! Init. Grille de proc. logiques. 88 'C', nprow, npcol) 89 90 call blacs_gridinfo(icntxt, & ! Mes coord. dans la grille 91 nprow, npcol, prow, pcol) ! de processeurs logiques. 92 93 call descinit(desc_a, n, n, & ! Init. descripteur de la 94 rb_size, cb_size, 0, 0, & ! matrice a 95 icntxt, n_max, info) 96 if (info.lt.0) stop 'descinit' 97 98 call descinit(desc_x, n, 1, & ! Init. descripteur du vecteur x 99 rb_size, 1, 0, 0, & 100 icntxt, n_max, info) & 101 if (info.lt.0) stop 'descinit' 102 103! Distribution de la matrice 'a' vers la matrice distribuee 'dist_a'. 104!-------------------------------------------------------------------- 105 call distribue(icntxt, a, n, dist_a, n_max, prow, pcol) 106 107! Distribution du vecteur 'b' sur le vecteur distribue 'dist_b'. 108!------------------------------------------------------------------ 109 if ((prow==0).and.(pcol==0)) then 110 dist_b(1:2) = b(1:2) ! copie de b(1:2) dans proc(0,0) 111 call sgesd2d(icntxt,2,1, & ! b(3:5) envoie au proc (1,0) 112 b(3),n,1,0) 113 dist_b(3) = b(5) ! copie de b(5) dans proc (0,0) 114 end if 115 if ((prow==1).and.(pcol==0)) then 116 call sgerv2d(icntxt, 2,1, dist_b(1), & ! reception b(3:5) dans 117 n_max, 0, 0) ! dist_b(1:2). 118 end if 119 120! Calculs 121!---------------------------------- 122 call PSGETRF(n, n, dist_a, 1, 1, desc_a, & 123 dist_piv, info) 124 if (info /= 0) print *,'Erreur dans la factorisation : ',info 125 126 call PSGETRS( 'N', n, 1, dist_a, 1, 1, desc_a, & 127 dist_piv, dist_b, 1, 1, desc_x, info) 128 if (info /= 0) print *,'Erreur dans la resolution : ',info 129 130! Recuperation des resultats 'dist_b' -> 'x' 131!---------------------------------------------- 132 if ((prow==1).and.(pcol==0)) then 133 call sgesd2d(icntxt,2,1,dist_b(1), & ! envoie de dist_b(1:2) 134 n_max, 0, 0) ! au proc (0,0) 135 end if 136 if ((prow==0).and.(pcol==0)) then 137 x(1:2) = dist_b(1:2) ! copie de dist_b(1:2) dans proc(0,0) 138 call sgerv2d(icntxt,2,1, & ! reception dist_b(1:2) dans x(3:4) 139 x(3), n, 1, 0) 140 x(5) = dist_b(3) ! copie de dist_b(3) dans proc (0,0) 141 end if 142 143! Impression des resultats 144!----------------------------- 145 if ((prow==0).and.(pcol==0)) then 146 write(6,*) ' resultats :' 147 write(6,'(1x,5F5.1)') (x(j),j=1,n) 148 end if 149 150end program 151 152!-------------------------------------------------------------------- 153! Distribution de la matrice 'x' vers la matrice distribuee 'dist_x'. 154!-------------------------------------------------------------------- 155subroutine distribue( icntxt, x, n, dist_x, n_max, prow, pcol) 156 implicit none 157 integer :: icntxt ! Contexte (grille de procs) 158 integer :: n ! Ordre de la matrice X 159 integer :: n_max ! Ordre de la matrice locale. 160 real,dimension(n,n):: x ! Matrice globale X. 161 real,dimension(n_max,n_max):: dist_x ! Matrice locale. 162 integer :: prow, pcol ! Coord. de mon proc. logique 163 164 if ((prow == 0) .and. (pcol == 0)) then ! Emission par le 165 ! processeur (0,0) 166 call sgesd2d( icntxt, 2,2, x(1,3), n, 0, 1) 167 call sgesd2d( icntxt, 2,2, x(3,1), n, 1, 0) 168 call sgesd2d( icntxt, 2,2, x(3,3), n, 1, 1) 169 call sgesd2d( icntxt, 2,1, x(3,5), n, 1, 0) 170 call sgesd2d( icntxt, 1,2, x(5,3), n, 0, 1) 171 end if 172 173 if ((prow == 0) .and. (pcol == 0)) then ! Copie des blocs 174 dist_x(1:2,1:2) = x(1:2,1:2) ! appartenant au proc(0,0) 175 dist_x(1:2,3) = x(1:2,5) 176 dist_x(3,1:2) = x(5,1:2) 177 dist_x(3,3) = x(5,5) 178 end if 179 ! Reception par les 180 if ((prow == 0) .and. (pcol == 1)) then ! autres processeurs. 181 call SGERV2D( icntxt, 2,2, dist_x(1,1), n_max, 0, 0) 182 call SGERV2D( icntxt, 1,2, dist_x(3,1), n_max, 0, 0) 183 end if 184 if ((prow == 1) .and. (pcol == 0)) then 185 call SGERV2D( icntxt, 2,2, dist_x(1,1), n_max, 0, 0) 186 call SGERV2D( icntxt, 2,1, dist_x(1,3), n_max, 0, 0) 187 end if 188 if ((prow == 1) .and. (pcol == 1)) then 189 call SGERV2D( icntxt, 2,2, dist_x(1,1), n_max, 0, 0) 190 end if 191 192end subroutine distribue 193