/* vide {rdo_sparse_matrix.h} */
#define rdo_sparse_matrix_C_COPYRIGHT "Copyright © 2008 Danillo Pereira e J. Stolfi, UNICAMP"
/* Last edited on 2008-07-22 17:41:45 by stolfi */ 

#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

#include <vec.h>
#include <affirm.h>
#include <filefmt.h>
#include <nget.h>
#include <fget.h>

#include "rdo_geometry.h"
#include "rdo_sparse_matrix.h"

MatrizEsparsa AlocaMatrizEsparsa(int nlin, int ncol, int nen)
  {
    MatrizEsparsa A;
    A.nlin = nlin;
    A.ncol = ncol;
    A.e = notnull(malloc(nen * sizeof(EntradaDeMatriz)), "falta memória");
    A.nen = nen;
    return A;
  }

void DesalocaMatrizEsparsa(MatrizEsparsa *A)
  {
    free(A->e); A->nen = 0;
  }

void ExpandeMatrizEsparsa(MatrizEsparsa *A, int pos)
  {
    vec_expand(&(A->nen), (void **)&(A->e), pos, sizeof(EntradaDeMatriz));
  }

void AjustaTamanhoDeMatrizEsparsa(MatrizEsparsa *A, int nen)
  {
    vec_trim(&(A->nen), (void **)&(A->e), nen, sizeof(EntradaDeMatriz));
  }

int AcrescentaEntradaAMatrizEsparsa(MatrizEsparsa *A, int pos, int lin, int col, double v)
  { if(v != 0.0)
      { ExpandeMatrizEsparsa(A, pos);
        A->e[pos] = (EntradaDeMatriz) { .lin = lin, .col = col, .valor = v };
        pos++;
      }
    return pos;
  }
  
int AcrescentaLinhaAMatrizEsparsa(MatrizEsparsa *A, int pos, int lin, double v[])
  {
    int nlin = A->nlin;
    int ncol = A->ncol;
    demand((lin >= 0) && (lin < nlin), "número de linha inválido");
    int j;
    for (j = 0; j < ncol; j++)
      { if(v[j] != 0.0)
	  { pos = AcrescentaEntradaAMatrizEsparsa(A, pos, lin, j, v[j]); }
      }
    return pos;
  }

void CopiaMatrizEsparsa(MatrizEsparsa *A, MatrizEsparsa *B)
  { AjustaTamanhoDeMatrizEsparsa(B, A->nen);
    int pos;
    for(pos = 0; pos < A->nen; pos++) { B->e[pos] = A->e[pos]; }
    B->nlin = A->nlin;
    B->ncol = A->ncol;
  }

void MultiplicaMatrizEsparsaPorVetor(MatrizEsparsa *A, double b[], double c[])
  { int nlin = A->nlin;
    int ncol = A->ncol;
    int i;
    for (i = 0; i < nlin; i++) { c[i] = 0.0; }
    int pos;
    for (pos = 0; pos < A->nen; pos++)
      { EntradaDeMatriz *eP = &(A->e[pos]);
        assert((eP->lin >= 0) && (eP->lin < nlin));
        assert((eP->col >= 0) && (eP->col < ncol));
        c[eP->lin] += eP->valor * b[eP->col];
      }
  }

void MultiplicaVetorPorMatrizEsparsa(double b[], MatrizEsparsa *A, double c[])
  { int nlin = A->nlin;
    int ncol = A->ncol;
    int j;
    for (j = 0; j < ncol; j++) { c[j] = 0.0; }
    int pos;
    for (pos = 0; pos < A->nen; pos++)
      { EntradaDeMatriz *eP = &(A->e[pos]); 
        assert((eP->lin >= 0) && (eP->lin < nlin));
        assert((eP->col >= 0) && (eP->col < ncol));
        c[eP->col] += eP->valor * b[eP->lin];
      }
  }

void MultiplicaMatrizesEsparsas(MatrizEsparsa *A, MatrizEsparsa *B, MatrizEsparsa *C)
  {
    demand(A->ncol == B->nlin, "matrizes de tamanhos incompatíveis");

    //aloca vetores de trabalho:
    double *av = malloc(A->ncol * sizeof(double)); //uma linha de {A}.
    double *cv = malloc(B->ncol * sizeof(double)); //uma linha de {C}.

    //multiplica cada linha da matriz {A} por {B}, insere em {C}: 
    int posA = 0;    //índice da primeira entrada da próxima linha de {A}.
    int nenC = 0;    //número de entradas já acrescentadas a {C}.
    while (posA < A->nen)
      { int linA = A->e[posA].lin; //índice da próxima linha na matriz {A} expandida.
        //coloca a próxima linha de {A} no vetor {av}:
        posA = ExtraiLinhaDeMatrizEsparsa(A, posA, linA, av);
       //multiplica essa linha de {A} por {B} colocando no vetor {cv}:
        MultiplicaVetorPorMatrizEsparsa(av, B, cv);
        //acrescenta a {C} os elementos não nulos de {cv}: 
        nenC = AcrescentaLinhaAMatrizEsparsa(C, nenC, linA, cv);
      }

    //corta {C} no tamanho justo:
    AjustaTamanhoDeMatrizEsparsa(C, nenC);
    //libera vetores de trabalho:
    free(av);
    free(cv);
  }

int ExtraiLinhaDeMatrizEsparsa(MatrizEsparsa *A, int pos, int lin, double v[])
  {
    int nlin = A->nlin;
    int ncol = A->ncol;
    demand((lin >= 0) && (lin  < nlin), "índice de linha inválido"); 
    
    //limpa o vetor {v}:
    int j;
    for (j = 0; j < ncol; j++) { v[j] = 0.0; }
    
    //extrai entradas dessa linha:
    while (pos < A->nen)
      { EntradaDeMatriz *eP = &(A->e[pos]); 
        assert((eP->lin >= 0) && (eP->lin < nlin));
        assert((eP->col >= 0) && (eP->col < ncol));
        if (eP->lin > lin) { break; }
        demand(eP->lin == lin, "entradas fora de ordem");
        v[eP->col] = eP->valor;
      }
    return pos;
  }
    
//identifica o formato atual da descrição da cena.
#define VERSAO_ARQ_MATRIZ "2008-07-21"

MatrizEsparsa LeMatrizEsparsa(FILE *rd)
  {
    MatrizEsparsa A;
    
    //lê cabecalho e eventuais comentários:
    filefmt_read_header(rd, "MatrizEsparsa", VERSAO_ARQ_MATRIZ);
    char *cmt = filefmt_read_comment(rd, '#');
    free(cmt); //por enquanto.

    //lê as dimensões:
    A.nlin = nget_uint(rd, "num_linhas", 10); fget_eol(rd);
    A.ncol = nget_uint(rd, "num_colunas", 10); fget_eol(rd);
    A.nen =  nget_uint(rd, "num_entradas", 10); fget_eol(rd);
    
    //lê as entradas:
    A.e = NOTNULL(malloc(A.nen * sizeof(EntradaDeMatriz)));
    int pos;
    for(pos = 0; pos < A.nen; pos++)
      {
        EntradaDeMatriz *eP = &(A.e[pos]);
        eP->lin = fget_uint(rd, 10);
        demand(eP->lin < A.nlin, "índice de linha inválido");
        eP->col = fget_uint(rd, 10);
        demand(eP->col < A.ncol, "índice de colun inválido");
        eP->valor = fget_double(rd);
        fget_eol(rd);    
      }
    
    //lê o rodapé:
    filefmt_read_footer(rd, "MatrizEsparsa");

    return A;
  }

void EscreveMatrizEsparsa(FILE *wr, MatrizEsparsa *A)
  {
    //escreve cabecalho:
    filefmt_write_header(wr, "MatrizEsparsa", VERSAO_ARQ_MATRIZ);
    //aqui deveria escrever comentários fornecidos pelo cliente.

    //escreve as dimensões:
    fprintf(wr, "num_linhas = %d\n", A->nlin);
    fprintf(wr, "num_colunas = %d\n", A->ncol);
    fprintf(wr, "num_entradas = %d\n", A->nen);
    
    //escreve as entradas:
    int pos;
    for(pos = 0; pos < A->nen; pos++)
      { EntradaDeMatriz *eP = &(A->e[pos]);
        fprintf(wr, "%5d %5d %23.16e\n", eP->lin, eP->col, eP->valor);
      }

    //escreve o rodapé e sincroniza arquivo:
    filefmt_write_footer(wr, "MatrizEsparsa");
    fflush(wr);
  }