#define RCSID "$Id: FMM_Operations.c,v 1.15 2006/02/26 00:42:53 geuzaine Exp $"
/*
* Copyright (C) 1997-2006 P. Dular, C. Geuzaine
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
* USA.
*
* Please report all bugs and problems to <getdp@geuz.org>.
*
* Contributor(s):
* Ruth Sabariego
*/
#include "GetDP.h"
#include "Numeric.h"
#include "Data_Passive.h"
#include "Data_DefineE.h"
#include "CurrentData.h"
#include "Cal_Value.h"
#include "Data_FMM.h"
/* ------------------------------------------------------------------------- */
/* FMM matrix vector product */
/* ------------------------------------------------------------------------- */
void FMM_MatVectorProd(double *x, double *y){
int NbrFMMEqu, iFMMEqu ;
struct FMMmat *FMMmat_P0, *FMMmat_P ;
GetDP_Begin("FMM_MatVectorProd");
NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
FMMmat_P = FMMmat_P0 + iFMMEqu ;
((void (*)(struct FMMmat*,double*,double*))FMMmat_P->FctProd)(FMMmat_P, x, y) ;
}
GetDP_End;
}
void FMM_NoRenumbering( ){
int NbrFMMEqu, iFMMEqu, iG, NbrGroupSrc, NbrGroupObs ;
struct FMMmat *FMMmat_P0, *FMMmat_P ;
int iDof, NbrDof, NumDof, *NumDof_A ;
int iEqu, NbrEqu, NumEqu, *NumEqu_A ;
int NbrHar, NumDofi, NumEqui, Flag_Iter = 0 ;
List_T *NumDof_L, *NumEqu_L ;
List_T *NumDofi_L, *NumEqui_L ;
GetDP_Begin("FMM_NoRenumbering");
NbrHar = Current.NbrHar;
NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
FMMmat_P = FMMmat_P0 + iFMMEqu ;
if(FMMmat_P->NumDofr == NULL)
FMMmat_P->NumDofr = FMMmat_P->NumDof ;
if(FMMmat_P->NumEqur == NULL)
FMMmat_P->NumEqur = FMMmat_P->NumEqu ;
if(NbrHar == 2){
NbrGroupSrc = List_Nbr( FMMmat_P->NumDof );
if(FMMmat_P->NumDofi == NULL)
FMMmat_P->NumDofi = List_Create(NbrGroupSrc, 1, sizeof(List_T*)) ;
else Flag_Iter = 1 ;
if(!Flag_Iter){
for (iG = 0 ; iG < NbrGroupSrc ; iG++){
List_Read(FMMmat_P->NumDof, iG, &NumDof_L ) ;
NbrDof = List_Nbr( NumDof_L);
NumDof_A = (int*)(NumDof_L->array) ;
NumDofi_L = List_Create(NbrDof, 1, sizeof(int)) ;
for (iDof = 0 ; iDof < NbrDof ; iDof++){
NumDof = NumDof_A[iDof] ;
NumDofi = NumDof+1 ;/* Imaginary part */
List_Add(NumDofi_L, &NumDofi) ;
}
List_Add(FMMmat_P->NumDofi, &NumDofi_L ) ;
}
if(FMMmat_P->NumEqu != FMMmat_P->NumDof){
NbrGroupObs = List_Nbr( FMMmat_P->NumEqu );
FMMmat_P->NumEqui = List_Create(NbrGroupObs, 1, sizeof(List_T*)) ;
for (iG = 0 ; iG < NbrGroupObs ; iG++){
List_Read(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
NbrEqu = List_Nbr( NumEqu_L);
NumEqu_A = (int*)(NumEqu_L->array) ;
NumEqui_L = List_Create(NbrEqu, 1, sizeof(int)) ;
for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
NumEqu = NumEqu_A[iEqu] ;
NumEqui = NumEqu+1 ;/* Imaginary part */
List_Add(NumEqui_L, &NumEqui) ;
}
List_Add(FMMmat_P->NumEqui, &NumEqui_L ) ;
}
}
else
FMMmat_P->NumEqui = FMMmat_P->NumDofi ;
}/* Flag_iter */
}/* NbrHar */
}/* iFMMEqu */
GetDP_End;
}
void FMM_Renumbering( int N, int *permr, int *permp){
int NbrFMMEqu, iFMMEqu, iG, NbrGroupSrc, NbrGroupObs ;
struct FMMmat *FMMmat_P0, *FMMmat_P ;
int iDof, NbrDof, NumDof, NumDofr, *NumDof_A ;
int iEqu, NbrEqu, NumEqu, NumEqur, *NumEqu_A ;
int NbrHar, NumDofri, NumEquri, Flag_Renum = 0 ;
List_T *NumDof_L, *NumDofr_L, *NumEqu_L , *NumEqur_L;
List_T *NumDofi_L, *NumEqui_L ;
GetDP_Begin("FMM_Renumbering");
NbrHar = Current.NbrHar;
NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
FMMmat_P = FMMmat_P0 + iFMMEqu ;
NbrGroupSrc = List_Nbr( FMMmat_P->NumDof );
if( FMMmat_P->NumDofr == NULL )
FMMmat_P->NumDofr = List_Create(NbrGroupSrc, 1, sizeof(List_T*)) ;
else Flag_Renum = 1 ;
if(NbrHar == 2 && FMMmat_P->NumDofi == NULL)
FMMmat_P->NumDofi = List_Create(NbrGroupSrc, 1, sizeof(List_T*)) ;
for (iG = 0 ; iG < NbrGroupSrc ; iG++){
List_Read(FMMmat_P->NumDof, iG, &NumDof_L ) ;
NbrDof = List_Nbr( NumDof_L);
NumDof_A = (int*)(NumDof_L->array) ;
if(!Flag_Renum)
NumDofr_L = List_Create(NbrDof, 1, sizeof(int)) ;
else
List_Read(FMMmat_P->NumDofr, iG, &NumDofr_L ) ;
if(NbrHar == 2){
if(!Flag_Renum)
NumDofi_L = List_Create(NbrDof, 1, sizeof(int)) ;
else
List_Read(FMMmat_P->NumDofi, iG, &NumDofi_L ) ;
}
for (iDof = 0 ; iDof < NbrDof ; iDof++){
NumDof = NumDof_A[iDof] - 1 ;
NumDofr = permp[(permr[NumDof]-1) + N] ;
if(!Flag_Renum)
List_Add(NumDofr_L, &NumDofr) ;
else
List_Write(NumDofr_L, iDof, &NumDofr) ;
if(NbrHar == 2){
NumDof++ ; /* Imaginary part */
NumDofri = permp[(permr[NumDof]-1) + N] ;
if(!Flag_Renum)
List_Add(NumDofi_L, &NumDofri) ;
else
List_Write(NumDofi_L, iDof, &NumDofri) ;
}
}
if(!Flag_Renum)
List_Add(FMMmat_P->NumDofr, &NumDofr_L ) ;
else
List_Write(FMMmat_P->NumDofr, iG, &NumDofr_L ) ;
if(NbrHar == 2){
if(!Flag_Renum)
List_Add(FMMmat_P->NumDofi, &NumDofi_L ) ;
else
List_Write(FMMmat_P->NumDofi, iG, &NumDofi_L ) ;
}
}
Flag_Renum = 0 ;
if(FMMmat_P->NumEqu != FMMmat_P->NumDof){
NbrGroupObs = List_Nbr( FMMmat_P->NumEqu );
if( FMMmat_P->NumEqur == NULL )
FMMmat_P->NumEqur = List_Create(NbrGroupSrc, 1, sizeof(List_T*)) ;
else Flag_Renum = 1 ;
if (NbrHar == 2 && FMMmat_P->NumEqui == NULL)
FMMmat_P->NumEqui = List_Create(NbrGroupObs, 1, sizeof(List_T*)) ;
for (iG = 0 ; iG < NbrGroupObs ; iG++){
List_Read(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
NbrEqu = List_Nbr( NumEqu_L);
NumEqu_A = (int*)(NumEqu_L->array) ;
if(!Flag_Renum)
NumEqur_L = List_Create(NbrEqu, 1, sizeof(int)) ;
else
List_Read(FMMmat_P->NumEqur, iG, &NumEqur_L ) ;
if(NbrHar == 2){
if(!Flag_Renum)
NumEqui_L = List_Create(NbrEqu, 1, sizeof(int)) ;
else
List_Read(FMMmat_P->NumEqui, iG, &NumEqui_L ) ;
}
for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
NumEqu = NumEqu_A[iEqu] - 1 ;
NumEqur = permp[(permr[NumEqu]-1) + N] ;
if(!Flag_Renum)
List_Add(NumEqur_L, &NumEqur) ;
else
List_Write(NumEqur_L, iEqu, &NumEqur) ;
if(NbrHar == 2){
NumEqu++ ; /* Imaginary part */
NumEquri = permp[(permr[NumEqu]-1) + N] ;
if(!Flag_Renum)
List_Add(NumEqui_L, &NumEquri) ;
else
List_Write(NumEqui_L, iEqu, &NumEquri) ;
}
}
if(!Flag_Renum)
List_Add(FMMmat_P->NumEqur, &NumEqur_L ) ;
else
List_Write(FMMmat_P->NumEqur, iG, &NumEqur_L ) ;
if(NbrHar == 2){
if(!Flag_Renum)
List_Add(FMMmat_P->NumEqui, &NumEqui_L ) ;
else
List_Write(FMMmat_P->NumEqui, iG, &NumEqui_L ) ;
}
}
}
else{
FMMmat_P->NumEqur = FMMmat_P->NumDofr ;
if (NbrHar == 2) FMMmat_P->NumEqui = FMMmat_P->NumDofi ;
}
Flag_Renum = 0;
}
GetDP_End;
}
void FMM_RenumberingOLD( int N, int *permr, int *permp){
int NbrFMMEqu, iFMMEqu, iG, NbrGroupSrc, NbrGroupObs ;
struct FMMmat *FMMmat_P0, *FMMmat_P ;
int iDof, NbrDof, NumDof, NumDofr, *NumDof_A ;
int iEqu, NbrEqu, NumEqu, NumEqur, *NumEqu_A ;
int NbrHar, NumDofri, NumEquri, Flag_Renum = 0 ;
List_T *NumDof_L, *NumEqu_L ;
List_T *NumDofi_L, *NumEqui_L ;
GetDP_Begin("FMM_Renumbering");
/* For using together with InverseRenumbering, the NumDof and NumEqu
are changed and for a second iteration they must be the original
ones. The Renumbering is done properly but not the
InverseRenumbering */
NbrHar = Current.NbrHar;
NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
FMMmat_P = FMMmat_P0 + iFMMEqu ;
NbrGroupSrc = List_Nbr( FMMmat_P->NumDof );
if(NbrHar == 2 && FMMmat_P->NumDofi == NULL)
FMMmat_P->NumDofi = List_Create(NbrGroupSrc, 1, sizeof(List_T*)) ;
else Flag_Renum = 1 ;
for (iG = 0 ; iG < NbrGroupSrc ; iG++){
List_Read(FMMmat_P->NumDof, iG, &NumDof_L) ;
NbrDof = List_Nbr( NumDof_L);
NumDof_A = (int*)(NumDof_L->array) ;
if(NbrHar == 2){
if(!Flag_Renum)
NumDofi_L = List_Create(NbrDof, 1, sizeof(int)) ;
else
List_Read(FMMmat_P->NumDofi, iG, &NumDofi_L ) ;
}
for (iDof = 0 ; iDof < NbrDof ; iDof++){
NumDof = NumDof_A[iDof] - 1 ;
NumDofr = permp[(permr[NumDof]-1) + N] ;
List_Write(NumDof_L, iDof, &NumDofr) ;
if(NbrHar == 2){
NumDof++ ; /* Imaginary part */
NumDofri = permp[(permr[NumDof]-1) + N] ;
if(!Flag_Renum)
List_Add(NumDofi_L, &NumDofri) ;
else
List_Write(NumDofi_L, iDof, &NumDofri) ;
}
}
List_Write(FMMmat_P->NumDof, iG, &NumDof_L ) ;
if(NbrHar == 2){
if(!Flag_Renum)
List_Add(FMMmat_P->NumDofi, &NumDofi_L ) ;
else
List_Write(FMMmat_P->NumDofi, iG, &NumDofi_L ) ;
}
}
Flag_Renum = 0 ;
if(FMMmat_P->NumEqu != FMMmat_P->NumDof){
NbrGroupObs = List_Nbr( FMMmat_P->NumEqu );
if (NbrHar == 2 && FMMmat_P->NumEqui == NULL)
FMMmat_P->NumEqui = List_Create(NbrGroupObs, 1, sizeof(List_T*)) ;
else Flag_Renum = 1 ;
for (iG = 0 ; iG < NbrGroupObs ; iG++){
List_Read(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
NbrEqu = List_Nbr( NumEqu_L);
NumEqu_A = (int*)(NumEqu_L->array) ;
if(NbrHar == 2){
if(!Flag_Renum)
NumEqui_L = List_Create(NbrEqu, 1, sizeof(int)) ;
else
List_Read(FMMmat_P->NumEqui, iG, &NumEqui_L ) ;
}
for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
NumEqu = NumEqu_A[iEqu] - 1 ;
NumEqur = permp[(permr[NumEqu]-1) + N] ;
List_Write(NumEqu_L, iEqu, &NumEqur) ;
if(NbrHar == 2){
NumEqu++ ; /* Imaginary part */
NumEquri = permp[(permr[NumEqu]-1) + N] ;
if(!Flag_Renum)
List_Add(NumEqui_L, &NumEquri) ;
else
List_Write(NumEqui_L, iEqu, &NumEquri) ;
}
}
List_Write(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
if(NbrHar == 2){
if(!Flag_Renum)
List_Add(FMMmat_P->NumEqui, &NumEqui_L ) ;
else
List_Write(FMMmat_P->NumEqui, iG, &NumEqui_L ) ;
}
}
}
else
if(NbrHar == 2) FMMmat_P->NumEqui = FMMmat_P->NumDofi ;
}
GetDP_End;
}
void FMM_InverseRenumbering( int *rpermr){
int NbrFMMEqu, iFMMEqu, iG, NbrGroupSrc, NbrGroupObs ;
struct FMMmat *FMMmat_P0, *FMMmat_P ;
int iDof, NbrDof, NumDof, NumDofr, *NumDof_A ;
int iEqu, NbrEqu, NumEqu, NumEqur, *NumEqu_A ;
List_T *NumDof_L, *NumEqu_L ;
GetDP_Begin("FMM_InverseRenumbering");
/* This function does not work properly! To revise... */
NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
FMMmat_P = FMMmat_P0 + iFMMEqu ;
NbrGroupSrc = List_Nbr( FMMmat_P->NumDof );
for (iG = 0 ; iG < NbrGroupSrc ; iG++){
List_Read(FMMmat_P->NumDof, iG, &NumDof_L ) ;
NbrDof = List_Nbr( NumDof_L);
NumDof_A = (int*)(NumDof_L->array) ;
for (iDof = 0 ; iDof < NbrDof ; iDof++){
NumDofr = NumDof_A[iDof] ;
NumDof = rpermr[NumDofr]-1 ;
List_Write(NumDof_L, iDof, &NumDof) ;
}
List_Write(FMMmat_P->NumDof, iG, &NumDof_L ) ;
}
if(FMMmat_P->NumEqu != FMMmat_P->NumDof){
NbrGroupObs = List_Nbr( FMMmat_P->NumEqu );
for (iG = 0 ; iG < NbrGroupObs ; iG++){
List_Read(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
NbrEqu = List_Nbr( NumEqu_L);
NumEqu_A = (int*)(NumEqu_L->array) ;
for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
NumEqur = NumEqu_A[iEqu] ;
NumEqu = rpermr[NumEqur]-1 ;
List_Write(NumEqu_L, iEqu, &NumEqu ) ;
}
List_Write(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
}
}
}
GetDP_End;
}
void FMM_Scaling(double *rowscal, double *colscal, int N){
int NbrFMMEqu, iFMMEqu, iG, NbrGroupSrc, NbrGroupObs;
int iR, iDir, NbrDir, iCom, NbrCom, Inc, TypeTimeDerivative, i ,Flag_scaled = 0 ;
struct FMMmat *FMMmat_P0, *FMMmat_P ;
int iDof, NbrDof, NumDof, *NumDof_A ;
int iEqu, NbrEqu, NumEqu, *NumEqu_A ;
double **Ag_M, **D_M, *Ag_V, *D_V, *x_prev = NULL ;
List_T *NumDof_L, *NumEqu_L ;
/* Aggregation and Disaggregation matrices can be either SCALAR or VECTOR */
GetDP_Begin("FMM_Scaling");
NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
NbrDir = Current.FMM.N ;
for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
FMMmat_P = FMMmat_P0 + iFMMEqu ;
TypeTimeDerivative = FMMmat_P->TypeTimeDerivative ;
if ( TypeTimeDerivative == DT_ && Current.NbrHar == 1 ){
#if defined(HAVE_SPARSKIT)
x_prev = ((Current.DofData->CurrentSolution-1)->x).V ;
#else
Msg(GERROR, "FMM_Scaling only implemented for Sparskit");
#endif
if (!Flag_scaled){
for (i = 0 ; i < N ; i++)
x_prev[i] /= colscal[i] ;
Flag_scaled = 1;
}
}
NbrCom = FMMmat_P->NbrCom ;
Inc = 2 * NbrCom ;
NbrGroupSrc = List_Nbr( FMMmat_P->NumDof );
for (iG = 0 ; iG < NbrGroupSrc ; iG++){
List_Read(FMMmat_P->A_L, iG, &Ag_M ) ;
List_Read(FMMmat_P->NumDof, iG, &NumDof_L ) ;
NbrDof = List_Nbr( NumDof_L);
NumDof_A = (int*)(NumDof_L->array) ;
for (iDof = 0 ; iDof < NbrDof ; iDof++){
NumDof = NumDof_A[iDof] - 1 ;
Ag_V = Ag_M[iDof] ;
for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=Inc){
for (iCom = 0 ; iCom < NbrCom ; iCom++){
Ag_V[iR+iCom] *= colscal[NumDof] ;
Ag_V[iR+NbrCom+iCom] *= colscal[NumDof] ;
}
}
}
}/* iG Source */
NbrGroupObs = List_Nbr( FMMmat_P->NumEqu );
for (iG = 0 ; iG < NbrGroupObs ; iG++){
List_Read(FMMmat_P->D_L, iG, &D_M ) ;
List_Read(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
NbrEqu = List_Nbr( NumEqu_L);
NumEqu_A = (int*)(NumEqu_L->array) ;
for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
NumEqu = NumEqu_A[iEqu] - 1 ;
D_V = D_M[iEqu] ;
for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=Inc){
for (iCom = 0 ; iCom < NbrCom ; iCom++){
D_V[iR+iCom] *= rowscal[NumEqu] ;
D_V[iR+NbrCom+iCom] *= rowscal[NumEqu] ;
}
}
}
}/* iG Observation */
}/* iFMMEqu */
GetDP_End;
}
void FMM_UnScaling(double *rowscal, double *colscal){
int NbrFMMEqu, iFMMEqu, iG, NbrGroupSrc, NbrGroupObs, iR, iDir, NbrDir, iCom, NbrCom, Inc ;
struct FMMmat *FMMmat_P0, *FMMmat_P ;
int iDof, NbrDof, NumDof, *NumDof_A ;
int iEqu, NbrEqu, NumEqu, *NumEqu_A ;
double **Ag_M, **D_M, *Ag_V, *D_V ;
List_T *NumDof_L, *NumEqu_L ;
/* Aggregation and Disaggregation matrices can be either SCALAR or VECTOR*/
GetDP_Begin("FMM_UnScaling");
NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
NbrDir = Current.FMM.N ;
for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
FMMmat_P = FMMmat_P0 + iFMMEqu ;
NbrCom = FMMmat_P->NbrCom ;
Inc = 2 * NbrCom ;
NbrGroupSrc = List_Nbr( FMMmat_P->NumDof );
for (iG = 0 ; iG < NbrGroupSrc ; iG++){
List_Read(FMMmat_P->A_L, iG, &Ag_M ) ;
List_Read(FMMmat_P->NumDof, iG, &NumDof_L ) ;
NbrDof = List_Nbr( NumDof_L);
NumDof_A = (int*)(NumDof_L->array) ;
for (iDof = 0 ; iDof < NbrDof ; iDof++){
NumDof = NumDof_A[iDof] - 1 ;
Ag_V = Ag_M[iDof] ;
for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=Inc){
for (iCom = 0 ; iCom < NbrCom ; iCom++){
Ag_V[iR+iCom] /= colscal[NumDof] ;
Ag_V[iR+NbrCom+iCom] /= colscal[NumDof] ;
}
}
}
}/* iG Source */
NbrGroupObs = List_Nbr( FMMmat_P->NumEqu );
for (iG = 0 ; iG < NbrGroupObs ; iG++){
List_Read(FMMmat_P->D_L, iG, &D_M ) ;
List_Read(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
NbrEqu = List_Nbr( NumEqu_L);
NumEqu_A = (int*)(NumEqu_L->array) ;
for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
NumEqu = NumEqu_A[iEqu] - 1 ;
D_V = D_M[iEqu] ;
for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=Inc){
for (iCom = 0 ; iCom < NbrCom ; iCom++){
D_V[iR+iCom] /= rowscal[NumEqu] ;
D_V[iR+NbrCom+iCom] /= rowscal[NumEqu] ;
}
}
}
}/* iG Observation */
}/* iFMMEqu */
GetDP_End;
}
void FMM_DTAMatrix(int N, double ***DTA){
int i, NbrFMMEqu, iFMMEqu ;
struct FMMmat *FMMmat_P0, *FMMmat_P ;
double **M, *x, *y ;
GetDP_Begin("FMM_DTAMatrix");
M =(double**)Malloc(N *sizeof(double*));
for(i = 0 ; i < N ; i++)
M[i] =(double*)Calloc(N,sizeof(double));
x = (double*)Calloc(N, sizeof(double));
NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
FMMmat_P = FMMmat_P0 + iFMMEqu ;
for (i = 0 ; i < N ; i++){
x[i] = 1. ;
y = M[i] ; /* row */
((void (*)(struct FMMmat*,double*,double*))FMMmat_P->FctProd)(FMMmat_P, x, y) ;
x[i] = 0. ;
}
}
*DTA = M ;
GetDP_End;
}
/* FMMProd_Laplace2D and FMMProd_Laplace3D admit only SCALAR Aggregation and Disaggregation,
the VECTOR case is not yet implemented and will be considered in separate routines */
void FMMProd_AllLaplace2DNonAdaptive(struct FMMmat *FMMmat_P, double *x, double *y ){
/* WARNING: Aggregation and Disaggregation matrices must be SCALAR */
int NbrGroupSrc, NbrDir ;
int NbrFG, *NumFG ;
int iDof, NbrDof, NumDof, *NumDof_A ;
int iEqu, NbrEqu, NumEqu, *NumEqu_A ;
int iG1, iG2, iDir, iDir2, iR, iHar, NbrHar ;
double **Ag_M, *Ag_V, AgJRe[NBR_MAX_DIR], AgJIm[NBR_MAX_DIR], ***T, *TG1G2, **Disag_M, *Disag_V ;
double Treal, Timag, TARe[NBR_MAX_DIR], TAIm[NBR_MAX_DIR] ;
List_T *FG_L, *NumDof_L, *NumEqu_L ;
GetDP_Begin("FMMProd_AllLaplace2DNonAdaptive");
T = FMMmat_P->T ;
NbrHar = Current.NbrHar ;
NbrDir = Current.FMM.NbrDir ;
NbrGroupSrc = List_Nbr( FMMmat_P->A_L );
for (iHar = 0 ; iHar < NbrHar ; iHar++){
for (iG1 = 0 ; iG1 < NbrGroupSrc ; iG1++){
List_Read(FMMmat_P->A_L, iG1, &Ag_M ) ;
for (iDir = 0; iDir < NbrDir ; iDir++)
AgJRe[iDir] = AgJIm[iDir] = 0.;
List_Read(FMMmat_P->FG_L, iG1, &FG_L) ;
NbrFG = List_Nbr(FG_L) ;
if(NbrFG){
if (iHar == 0)
List_Read(FMMmat_P->NumDof, iG1, &NumDof_L ) ;
else
List_Read(FMMmat_P->NumDofi, iG1, &NumDof_L ) ;
NbrDof = List_Nbr(NumDof_L);
NumDof_A = (int*)(NumDof_L->array) ;
for (iDof = 0 ; iDof < NbrDof ; iDof++){
NumDof = NumDof_A[iDof] - 1 ;
Ag_V = Ag_M[iDof] ;
for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=2){
AgJRe[iDir] += Ag_V[iR ] * x[NumDof] ;
AgJIm[iDir] += Ag_V[iR+1] * x[NumDof] ;
}
} /* iDof NbrDof */
NumFG = (int*)(FG_L->array) ;
for(iG2=0 ; iG2 < NbrFG ; iG2++){ /* Destination */
List_Read(FMMmat_P->D_L, NumFG[iG2], &Disag_M ) ;
TG1G2 = T[iG1][iG2] ;
for (iDir = 0, iR = 0 ; iDir < NbrDir ; iDir++){ /* Translation is not diagonal! */
TARe[iDir] = TAIm[iDir] = 0. ;
for (iDir2 = 0 ; iDir2 < NbrDir ; iDir2++, iR+=2){
Treal = TG1G2[iR ] ;
Timag = TG1G2[iR+1] ;
TARe[iDir] += Treal * AgJRe[iDir2] - Timag * AgJIm[iDir2] ;
TAIm[iDir] += Treal * AgJIm[iDir2] + Timag * AgJRe[iDir2] ;
}/* iDir2 */
}/* iDir1 */
if(iHar == 0)
List_Read(FMMmat_P->NumEqu, NumFG[iG2], &NumEqu_L ) ;
else
List_Read(FMMmat_P->NumEqui, NumFG[iG2], &NumEqu_L ) ;
NbrEqu = List_Nbr( NumEqu_L ) ;
NumEqu_A = (int*)( NumEqu_L->array) ;
for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
NumEqu = NumEqu_A[iEqu] - 1 ;
Disag_V = Disag_M[iEqu] ;
for (iDir = 0, iR = 0 ; iDir < NbrDir ; iDir++, iR+=2)
y[NumEqu] += Disag_V[iR] * TARe[iDir] - Disag_V[iR+1] * TAIm[iDir] ;
}/* iEqu NbrEqu*/
}/* iG2 Destination*/
}/* if NbrFG */
}/* iG1 Origin */
}/* iHar */
GetDP_End ;
}
void FMMProd_AllLaplace2D(struct FMMmat *FMMmat_P, double *x, double *y ){
/* WARNING: Aggregation and Disaggregation matrices must be SCALAR */
int NbrGroupSrc, NbrDir ;
int NbrFG, *NumFG, *Nd_A ;
int iDof, NbrDof, NumDof, NumDofprev, *NumDof_A, *NumDofprev_A = NULL;
int iEqu, NbrEqu, NumEqu, *NumEqu_A ;
int iG1, iG2, iDir, iDir2, iR, iHar, NbrHar,TypeTimeDerivative ;
double **Ag_M, *Ag_V, AgJRe[NBR_MAX_DIR], AgJIm[NBR_MAX_DIR], ***T, *TG1G2, **Disag_M, *Disag_V ;
double Treal, Timag, TARe[NBR_MAX_DIR], TAIm[NBR_MAX_DIR], w, *x_prev = NULL;
List_T *FG_L, *Nd_L, *NumDof_L, *NumDofprev_L, *NumEqu_L ;
GetDP_Begin("FMMProd_AllLaplace2D");
T = FMMmat_P->T ;
NbrHar = Current.NbrHar ;
TypeTimeDerivative = FMMmat_P->TypeTimeDerivative ;
w = FMMmat_P->Pulsation ;
if (TypeTimeDerivative == DT_ && NbrHar == 1 )
#if defined(HAVE_SPARSKIT)
x_prev = ((Current.DofData->CurrentSolution-1)->x).V ;
#else
Msg(GERROR, "FMM_Scaling only implemented for Sparskit");
#endif
NbrGroupSrc = List_Nbr( FMMmat_P->A_L );
for (iHar = 0 ; iHar < NbrHar ; iHar++){
for (iG1 = 0 ; iG1 < NbrGroupSrc ; iG1++){
List_Read(FMMmat_P->A_L, iG1, &Ag_M ) ;
NbrDir = Current.FMM.NbrDir ;
for (iDir = 0; iDir < NbrDir ; iDir++)
AgJRe[iDir] = AgJIm[iDir] = 0.;
List_Read(FMMmat_P->FG_L, iG1, &FG_L) ;
NbrFG = List_Nbr(FG_L) ;
if(NbrFG){
List_Read(FMMmat_P->Nd_L, iG1, &Nd_L) ;
Nd_A = (int*)(Nd_L->array) ;
if (iHar == 0)
List_Read(FMMmat_P->NumDofr, iG1, &NumDof_L ) ;
else
List_Read(FMMmat_P->NumDofi, iG1, &NumDof_L ) ;
if(TypeTimeDerivative == DT_ && NbrHar ==1 ){
List_Read(FMMmat_P->NumDof, iG1, &NumDofprev_L ) ;
NumDofprev_A = (int*)(NumDofprev_L->array) ;
}
NbrDof = List_Nbr(NumDof_L);
NumDof_A = (int*)(NumDof_L->array) ;
for (iDof = 0 ; iDof < NbrDof ; iDof++){
NumDof = NumDof_A[iDof] - 1 ;
Ag_V = Ag_M[iDof] ;
if(TypeTimeDerivative == NODT_ || ( TypeTimeDerivative == DT_ && NbrHar == 2) ){
for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=2){
AgJRe[iDir] += Ag_V[iR ] * x[NumDof] ;
AgJIm[iDir] += Ag_V[iR+1] * x[NumDof] ;
}
}
if (TypeTimeDerivative == DT_ && NbrHar ==1 ){
NumDofprev = NumDofprev_A[iDof] - 1 ;
for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=2){
AgJRe[iDir] += Ag_V[iR ] * ( x[NumDof] - x_prev[NumDofprev] ) ;
AgJIm[iDir] += Ag_V[iR+1] * ( x[NumDof] - x_prev[NumDofprev] ) ;
}
}
} /* iDof NbrDof */
NumFG = (int*)(FG_L->array) ;
for(iG2=0 ; iG2 < NbrFG ; iG2++){ /* Destination */
NbrDir = Nd_A[iG2] ;
List_Read(FMMmat_P->D_L, NumFG[iG2], &Disag_M ) ;
TG1G2 = T[iG1][iG2] ;
for (iDir = 0, iR = 0 ; iDir < NbrDir ; iDir++){ /* Translation is not diagonal! */
TARe[iDir] = TAIm[iDir] = 0. ;
for (iDir2 = 0 ; iDir2 < NbrDir ; iDir2++, iR+=2){
Treal = TG1G2[iR ] ;
Timag = TG1G2[iR+1] ;
TARe[iDir] += Treal * AgJRe[iDir2] - Timag * AgJIm[iDir2] ;
TAIm[iDir] += Treal * AgJIm[iDir2] + Timag * AgJRe[iDir2] ;
}/* iDir2 */
}/* iDir1 */
if(iHar == 0)
List_Read(FMMmat_P->NumEqur, NumFG[iG2], &NumEqu_L ) ;
else
List_Read(FMMmat_P->NumEqui, NumFG[iG2], &NumEqu_L ) ;
NbrEqu = List_Nbr( NumEqu_L ) ;
NumEqu_A = (int*)( NumEqu_L->array) ;
for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
NumEqu = NumEqu_A[iEqu] - 1 ;
Disag_V = Disag_M[iEqu] ;
for (iDir = 0, iR = 0 ; iDir < NbrDir ; iDir++, iR+=2){
y[NumEqu] += Disag_V[iR] * TARe[iDir] - Disag_V[iR+1] * TAIm[iDir] ;
}
}/* iEqu NbrEqu*/
}/* iG2 Destination*/
}/* if NbrFG */
}/* iG1 Origin */
}/* iHar */
GetDP_End ;
}
void FMMProd_AllLaplace3D(struct FMMmat *FMMmat_P, double *x, double *y ){
/* Aggregation and Disaggregation matrices SCALAR */
int NbrGroupSrc, NbrDir ;
int NbrFG, *NumFG, *Nd_A ;
int iDof, NbrDof, NumDof, NumDofprev, *NumDof_A, *NumDofprev_A = NULL;
int iEqu, NbrEqu, NumEqu, *NumEqu_A ;
int iG1, iG2, iDir, iDir2, iR, iHar, NbrHar, l, m, j, k ;
int TypeTimeDerivative ;
double **Ag_M, *Ag_V, AgJRe[NBR_MAX_DIR], AgJIm[NBR_MAX_DIR], ***T, *TG1G2, **Disag_M, *Disag_V ;
double Treal, Timag, TARe[NBR_MAX_DIR], TAIm[NBR_MAX_DIR], w, *x_prev = NULL ;
List_T *FG_L, *Nd_L, *NumDof_L, *NumDofprev_L, *NumEqu_L ;
GetDP_Begin("FMMProd_AllLaplace3D");
T = FMMmat_P->T ;
NbrHar = Current.NbrHar ;
TypeTimeDerivative = FMMmat_P->TypeTimeDerivative ;
w = FMMmat_P->Pulsation ;
if (TypeTimeDerivative == DT_ && NbrHar == 1 )
#if defined(HAVE_SPARSKIT)
x_prev = ((Current.DofData->CurrentSolution-1)->x).V ;
#else
Msg(GERROR, "FMM_Scaling only implemented for Sparskit");
#endif
NbrGroupSrc = List_Nbr( FMMmat_P->A_L );
for (iHar = 0 ; iHar < NbrHar ; iHar++){
for (iG1 = 0 ; iG1 < NbrGroupSrc ; iG1++){
List_Read(FMMmat_P->A_L, iG1, &Ag_M ) ;
NbrDir = (Current.FMM.NbrDir-1)*(Current.FMM.NbrDir+1) ;
for (iDir = 0; iDir < NbrDir ; iDir++)
AgJRe[iDir] = AgJIm[iDir] = 0.;
List_Read(FMMmat_P->FG_L, iG1, &FG_L) ;
NbrFG = List_Nbr(FG_L) ;
if(NbrFG){
List_Read(FMMmat_P->Nd_L, iG1, &Nd_L) ;
Nd_A = (int*)(Nd_L->array) ;
if (iHar == 0)
List_Read(FMMmat_P->NumDofr, iG1, &NumDof_L ) ;
else
List_Read(FMMmat_P->NumDofi, iG1, &NumDof_L ) ;
if(TypeTimeDerivative == DT_ && NbrHar ==1 ){
List_Read(FMMmat_P->NumDof, iG1, &NumDofprev_L ) ;
NumDofprev_A = (int*)(NumDofprev_L->array) ;
}
NbrDof = List_Nbr(NumDof_L);
NumDof_A = (int*)(NumDof_L->array) ;
for (iDof = 0 ; iDof < NbrDof ; iDof++){
NumDof = NumDof_A[iDof] - 1 ;
Ag_V = Ag_M[iDof] ;
if(TypeTimeDerivative == NODT_ || ( TypeTimeDerivative == DT_ && NbrHar == 2) )
for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=2){
AgJRe[iDir] += Ag_V[iR ] * x[NumDof] ;
AgJIm[iDir] += Ag_V[iR+1] * x[NumDof] ;
}
if (TypeTimeDerivative == DT_ && NbrHar ==1 ){
NumDofprev = NumDofprev_A[iDof] - 1 ;
for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=2){
AgJRe[iDir] += Ag_V[iR ] * ( x[NumDof] - x_prev[NumDofprev] ) ;
AgJIm[iDir] += Ag_V[iR+1] * ( x[NumDof] - x_prev[NumDofprev] ) ;
}
}
} /*iDof NbrDof*/
NumFG = (int*)(FG_L->array) ;
for(iG2=0 ; iG2 < NbrFG ; iG2++){ /* Destination */
NbrDir = (Nd_A[iG2]-1)*(Nd_A[iG2]+1) ;
List_Read(FMMmat_P->D_L, NumFG[iG2], &Disag_M ) ;
TG1G2 = T[iG1][iG2] ;
l = 0 ; m = 0 ;
for (iDir = 0; iDir < NbrDir ; iDir++){
j = 0 ; k = 0 ;
TARe[iDir] = TAIm[iDir] = 0 ;
for (iDir2 = 0 ; iDir2 < NbrDir ; iDir2++){
iR = 2*((l+j)*(l+j+1)+m+k) ;
Treal = TG1G2[iR ] ;
Timag = TG1G2[iR+1] ;
TARe[iDir] += Treal * AgJRe[iDir2] - Timag * AgJIm[iDir2] ;
TAIm[iDir] += Treal * AgJIm[iDir2] + Timag * AgJRe[iDir2] ;
if (j==k) {j++;k=-j;} else k++;
}/* iDir2 NbrDir */
if (l==m) {l++;m=-l;} else m++;
}/* iDir NbrDir */
if ( (TypeTimeDerivative == NODT_) || (TypeTimeDerivative == DT_ && NbrHar == 1) ){
if(iHar == 0)
List_Read(FMMmat_P->NumEqur, NumFG[iG2], &NumEqu_L ) ;
else
List_Read(FMMmat_P->NumEqui, NumFG[iG2], &NumEqu_L ) ;
}
if (TypeTimeDerivative == DT_ && NbrHar == 2){
if (iHar == 0)
List_Read(FMMmat_P->NumEqui, NumFG[iG2], &NumEqu_L ) ;
else
List_Read(FMMmat_P->NumEqur, NumFG[iG2], &NumEqu_L ) ;
}
NbrEqu = List_Nbr( NumEqu_L ) ;
NumEqu_A = (int*)(NumEqu_L->array) ;
for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
NumEqu = NumEqu_A[iEqu] - 1 ;
Disag_V = Disag_M[iEqu] ;
for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=2)
if ((TypeTimeDerivative == NODT_)|| (TypeTimeDerivative == DT_ && NbrHar == 1))
y[NumEqu] += Disag_V[iR] * TARe[iDir] - Disag_V[iR+1] * TAIm[iDir] ;
else
if (iHar ==0)
y[NumEqu] += w * (Disag_V[iR] * TARe[iDir] - Disag_V[iR+1] * TAIm[iDir]) ;
else
y[NumEqu] -= w * (Disag_V[iR] * TARe[iDir] - Disag_V[iR+1] * TAIm[iDir]) ;
}/* iEqu NbrEqu */
}/* iG2 Destination */
}/* if NbrFG */
}/* iG1 Origin */
}/* iHar */
GetDP_End ;
}
void FMMProd_AllLaplace3D_VECTOR(struct FMMmat *FMMmat_P, double *x, double *y ){
/* Aggregation and Disaggregation matrices can be either SCALAR or VECTOR*/
int NbrGroupSrc, NbrDir ;
int NbrFG, *NumFG, *Nd_A ;
int iDof, NbrDof, NumDof, *NumDof_A ;
int iEqu, NbrEqu, NumEqu, NumEquY, *NumEqu_A, *NumEquY_A ;
int iG1, iG2, iDir, iDir2, iR, iHar, NbrHar, l, m, j, k, iCom, NbrCom,Inc,iA,iTA ;
int TypeTimeDerivative ;
double **Ag_M, *Ag_V, AgJRe[NBR_MAX_DIR], AgJIm[NBR_MAX_DIR], ***T, *TG1G2, **Disag_M, *Disag_V ;
double Treal, Timag, TARe[NBR_MAX_DIR], TAIm[NBR_MAX_DIR], w ;
List_T *FG_L, *Nd_L, *NumDof_L, *NumEqu_L, *NumEquY_L ;
GetDP_Begin("FMMProd_AllLaplace3D");
T = FMMmat_P->T ;
TypeTimeDerivative = FMMmat_P->TypeTimeDerivative ;
w = FMMmat_P->Pulsation ;
NbrHar = Current.NbrHar ;
NbrCom = FMMmat_P->NbrCom ;
Inc = 2*NbrCom;
NbrGroupSrc = List_Nbr( FMMmat_P->A_L );
for (iHar = 0 ; iHar < NbrHar ; iHar++){
for (iG1 = 0 ; iG1 < NbrGroupSrc ; iG1++){
List_Read(FMMmat_P->A_L, iG1, &Ag_M ) ;
NbrDir = (Current.FMM.NbrDir-1)*(Current.FMM.NbrDir+1) ;
for (iDir = 0; iDir < NbrCom* NbrDir ; iDir++)
AgJRe[iDir] = AgJIm[iDir] = 0.;
List_Read(FMMmat_P->FG_L, iG1, &FG_L) ;
NbrFG = List_Nbr(FG_L) ;
if(NbrFG){
List_Read(FMMmat_P->Nd_L, iG1, &Nd_L) ;
Nd_A = (int*)(Nd_L->array) ;
if (iHar == 0)
List_Read(FMMmat_P->NumDof, iG1, &NumDof_L ) ;
else
List_Read(FMMmat_P->NumDofi, iG1, &NumDof_L ) ;
NbrDof = List_Nbr(NumDof_L);
NumDof_A = (int*)(NumDof_L->array) ;
for (iDof = 0 ; iDof < NbrDof ; iDof++){
NumDof = NumDof_A[iDof] - 1 ;
Ag_V = Ag_M[iDof] ;
for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=Inc){
for (iCom = 0 ; iCom < NbrCom ; iCom++){
AgJRe[iDir*NbrCom+iCom] += Ag_V[iR+iCom] * x[NumDof] ;
AgJIm[iDir*NbrCom+iCom] += Ag_V[iR+NbrCom+iCom] * x[NumDof] ;
}
}
} /*iDof NbrDof*/
NumFG = (int*)(FG_L->array) ;
for(iG2=0 ; iG2 < NbrFG ; iG2++){ /* Destination */
NbrDir = (Nd_A[iG2]-1)*(Nd_A[iG2]+1) ;
List_Read(FMMmat_P->D_L, NumFG[iG2], &Disag_M ) ;
TG1G2 = T[iG1][iG2] ;
l = 0 ; m = 0 ;
for (iDir = 0; iDir < NbrDir ; iDir++){
j = 0 ; k = 0 ;
TARe[iDir] = TAIm[iDir] = 0 ;
for (iDir2 = 0 ; iDir2 < NbrDir ; iDir2++){
iR = 2*((l+j)*(l+j+1)+m+k) ;
Treal = TG1G2[iR ] ;
Timag = TG1G2[iR+1] ;
for (iCom = 0 ; iCom < NbrCom ; iCom++){
iTA = iDir*NbrCom+iCom ;
iA = iDir2*NbrCom+iCom ;
TARe[iTA] += Treal * AgJRe[iA] - Timag * AgJIm[iA] ;
TAIm[iTA] += Treal * AgJIm[iA] + Timag * AgJRe[iA] ;
}
if (j==k) {j++;k=-j;} else k++;
}/* iDir2 NbrDir */
if (l==m) {l++;m=-l;} else m++;
}/* iDir NbrDir */
if(iHar == 0)
List_Read(FMMmat_P->NumEqu, NumFG[iG2], &NumEqu_L ) ;
else
List_Read(FMMmat_P->NumEqui, NumFG[iG2], &NumEqu_L ) ;
if (TypeTimeDerivative == NODT_)
NumEquY_L = NumEqu_L ;
else
if (iHar == 0)
List_Read(FMMmat_P->NumEqui, NumFG[iG2], &NumEquY_L ) ;
else
List_Read(FMMmat_P->NumEqu, NumFG[iG2], &NumEquY_L ) ;
NbrEqu = List_Nbr( NumEqu_L ) ;
NumEqu_A = (int*)(NumEqu_L->array) ;
NumEquY_A = (int*)(NumEquY_L->array) ;
for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
NumEqu = NumEqu_A[iEqu] - 1 ;
NumEquY = NumEquY_A[iEqu] - 1 ;
Disag_V = Disag_M[iEqu] ;
for (iDir = 0, iR = 0; iDir < NbrDir ; iDir++, iR+=Inc)
for (iCom = 0 ; iCom < NbrCom ; iCom++){
iTA = iDir*NbrCom+iCom ;
if (TypeTimeDerivative == NODT_)
y[NumEqu] += Disag_V[iR+iCom] * TARe[iTA] - Disag_V[iR+NbrCom+iCom] * TAIm[iTA] ;
else
if (iHar ==0)
y[NumEquY] -= w * (Disag_V[iR+iCom] * TARe[iTA] - Disag_V[iR+NbrCom+iCom] * TAIm[iTA]) ;
else
y[NumEquY] += w * (Disag_V[iR+iCom] * TARe[iTA] - Disag_V[iR+NbrCom+iCom] * TAIm[iTA]) ;
}
}/* iEqu NbrEqu */
}/* iG2 Destination */
}/* if NbrFG */
}/* iG1 Origin */
}/* iHar */
GetDP_End ;
}
void FMMProd_AllHelmholtz(struct FMMmat *FMMmat_P, double *x, double *y ){
int NbrGroupSrc, NbrDir;
int NbrFG, *NumFG ;
int iDof, NbrDof, NumDof, *NumDof_A ;
int iEqu, NbrEqu, NumEqu, NumEqui, *NumEqu_A, *NumEqui_A ;
int iG1, iG2, iDir, iHar, NbrHar, iR, iCom, NbrCom,iA,iRD ;
double **Ag_M, *Ag_V, AgJRe[3*NBR_MAX_DIR], AgJIm[3*NBR_MAX_DIR], ***T, *TG1G2, **Disag_M, *Disag_V ;
double TrAr, TrAi, TiAi, TiAr ;
List_T *FG_L, *NumDof_L, *NumEqu_L,*NumEqui_L ;
GetDP_Begin("FMMProd_AllHelmholtz");
NbrDir = Current.FMM.NbrDir ;
NbrHar = Current.NbrHar ;
NbrCom = FMMmat_P->NbrCom ;
T = FMMmat_P->T ;
NbrGroupSrc = List_Nbr( FMMmat_P->A_L );
for (iHar = 0 ; iHar < NbrHar ; iHar++){
for (iG1 = 0 ; iG1 < NbrGroupSrc ; iG1++){
List_Read(FMMmat_P->A_L, iG1, &Ag_M ) ;
for (iDir = 0; iDir < NbrDir ; iDir++)
for (iCom = 0 ; iCom < NbrCom ; iCom++)
AgJRe[iDir*NbrCom+iCom] = AgJIm[iDir*NbrCom+iCom] = 0 ;
List_Read(FMMmat_P->FG_L, iG1, &FG_L) ;
NbrFG = List_Nbr(FG_L) ;
if(NbrFG){
if(iHar == 0)
List_Read(FMMmat_P->NumDofr,iG1, &NumDof_L) ;
else
List_Read(FMMmat_P->NumDofi,iG1, &NumDof_L) ;
NbrDof = List_Nbr(NumDof_L);
NumDof_A = (int*)(NumDof_L->array) ;
for (iDof = 0 ; iDof < NbrDof ; iDof++){
NumDof = NumDof_A[iDof] - 1 ;
Ag_V = Ag_M[iDof] ;
for (iDir = 0, iR = 0; iDir < NbrDir; iDir++, iR+=2*NbrCom){
for (iCom = 0 ; iCom < NbrCom ; iCom++){
if(iHar == 0){
AgJRe[iDir*NbrCom+iCom] += Ag_V[iR+iCom ] * x[NumDof];
AgJIm[iDir*NbrCom+iCom] += Ag_V[iR+NbrCom+iCom] * x[NumDof];
}
else{
AgJRe[iDir*NbrCom+iCom] -= Ag_V[iR+NbrCom+iCom] * x[NumDof];
AgJIm[iDir*NbrCom+iCom] += Ag_V[iR+iCom] * x[NumDof];
}
}
}
}
NumFG = (int*)(FG_L->array) ;
for(iG2=0 ; iG2 < NbrFG ; iG2++){ /* Destination */
List_Read(FMMmat_P->D_L, NumFG[iG2], &Disag_M ) ;
TG1G2 = T[iG1][iG2] ;
List_Read(FMMmat_P->NumEqur, NumFG[iG2], &NumEqu_L) ;
List_Read(FMMmat_P->NumEqui, NumFG[iG2], &NumEqui_L) ;
NbrEqu = List_Nbr( NumEqu_L ) ;
NumEqu_A = (int*)( NumEqu_L->array ) ;
NumEqui_A = (int*)( NumEqui_L->array ) ;
for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
NumEqu = NumEqu_A[iEqu] - 1 ;
NumEqui = NumEqui_A[iEqu] - 1 ;
Disag_V = Disag_M[iEqu] ;
for (iDir = 0, iR = 0, iRD = 0 ; iDir < NbrDir ; iDir++, iR +=2, iRD+=2*NbrCom){
for (iCom = 0 ; iCom < NbrCom ; iCom++){
iA = iDir* NbrCom +iCom ;
TrAr = TG1G2[iR] * AgJRe[iA] ;
TrAi = TG1G2[iR] * AgJIm[iA] ;
TiAi = TG1G2[iR+1] * AgJIm[iA] ;
TiAr = TG1G2[iR+1] * AgJRe[iA] ;
y[NumEqu] += Disag_V[iRD+iCom] * ( TrAr - TiAi ) - Disag_V[iRD+NbrCom+iCom] * ( TrAi + TiAr ) ;
y[NumEqui] += Disag_V[iRD+iCom] * ( TrAi + TiAr ) + Disag_V[iRD+NbrCom+iCom] * ( TrAr - TiAi ) ;
}
}
}/* for iEqu NbrEqu*/
}/* for iG2 Destination*/
}/* if NbrFG */
}/* for iG1 Origin */
}/* for iHar */
GetDP_End ;
}
syntax highlighted by Code2HTML, v. 0.9.1